Method and apparatus for modeling injection of a fluid in a mold cavity

ABSTRACT

The invention relates to a method and apparatus for analyzing fluid flow while considering heat transfer effects and, in particular, a phase change from a molten state to a solid state. In particular, the method and apparatus may be applied to the analysis of an injection molding process for producing a molded polymer component from a thermoplastic or a thermosetting polymer. In one embodiment, the method may be used to determine pressure required to fill a mold cavity and pressure gradients introduced during filling and packing of the cavity of an injection mold. The results of these analyses may be used to determine the number and location of gates, to determine the best material for the component, and to optimize the process conditions used in the molding process.

TECHNICAL FIELD

The invention relates to the field of three dimensional modeling offluid flow in a cavity and, more specifically in one embodiment, to themodeling of an injection molding process for producing molded polymercomponents.

BACKGROUND

The use of injection molded plastic components has dramaticallyincreased in many industries in recent years. Manufacturers ofelectronic equipment, consumer goods, medical equipment, and automotiveparts are producing more and more of their products and components usedin their products out of plastics than ever before. At the same time,competitive pressures are driving manufacturers in the plasticsinjection molding industry to find new methods to optimize the designsin order to better match the designs to the production process. When theneed for component or mold configuration modifications are discoveredlate in the design development process, the delay and associated coststo implement the necessary changes rise rapidly. Companies that want toensure that their components are producible and will perform optimallyhave begun to use computer aided engineering techniques to simulate ormodel the complex flows in an injection mold, in order to understandbetter the manufacturing process and integrate this knowledge intocomponent design, early in the design phase.

There are a number of factors that should be considered when designingan injection mold and the component which is to be produced therein.Parameters such as the overall component geometry, minimum and maximumwall thicknesses, number and location of gates in the mold through whichthe liquid polymer is injected, number and location of vents in the moldthrough which gas in the cavity escapes, polymer composition andproperties, and shrinkage allowances, are a few. Due to the closelyinterrelated relationship, component and mold design cannot reliably bebased purely on form and function of the end component, but should alsoconsider the effects of the manufacturing process.

Computer aided engineering simulation can be used advantageously toprovide design and manufacturing engineers with visual and numericalfeedback as to what is likely to happen inside the mold cavity duringthe injection molding process, allowing them to better understand andpredict the behavior of contemplated component designs so that thetraditional, costly trial and error approach to manufacturing can beeliminated substantially. The use of computer aided engineeringsimulation facilitates optimizing component designs, mold designs, andmanufacturing processing parameters during the design phase, wherenecessary changes can be implemented readily, with the least cost andimpact on schedule.

A basic discussion of the injection molding process and the challengesassociated with producing high yield, quality injection moldedcomponents is addressed in the primer “Moldflow Design Principles:Quality and Productivity by Design,” distributed by Moldflow Pty. Ltd.,Kilsyth, Victoria, Australia, the assignee of the instant patentapplication, the disclosure of which is herein incorporated by referencein its entirety.

Briefly, the injection molding process is a complex, two step process.In the first step, referred to as the filling phase, polymer material isforced under pressure into the mold cavity until the cavity isvolumetrically filled. Thereafter, in the second step, referred to asthe packing phase, pressure is maintained on the polymer to permitfurther flow of polymer into the cavity to compensate for shrinkage asthe material solidifies and contracts. When the component issufficiently solid, the component may be ejected from the mold. Boththermoplastic and thermosetting polymers can be injection molded.

When molding thermoplastics, the temperature of the mold at the cavitysurface or wall is maintained at a temperature below the meltingtemperature of the material to be injected. As the material flows intothe cavity, the liquid material forms a solidified layer along thecavity wall. This layer may be referred to as the frozen layer and,depending on the processing conditions and the material used, thethickness thereof may vary during filling. The thickness of the frozenlayer is important, because the frozen layer reduces the effectivechannel width for flow in the cavity and, due to the thermo-rheologicalcharacteristics of thermoplastics, typically effects the viscosity ofthe material flowing thereby.

Early analytical simulation techniques relied upon two dimensionalfinite element models, which were found to be beneficial for simulatinginjection molding of relatively simple, thin-walled components. Moreadvanced simulation techniques are discussed, for example, inInternational Patent Application No. PCT/AU98/00130, assigned to theassignee of the instant invention, the disclosure of which isincorporated herein by reference in its entirety. In thick or complexcomponents, however, where molten plastic can flow in all directions,traditional thin-wall analytic assumptions that rely on planar regionsof specified thickness are not capable typically of predicting this typeof flow. To achieve high accuracy and predictability, a full, threedimensional simulation is desirable to establish, for example, whereweld lines will form, air traps will occur, and flow will lead or lag.

In order to analyze in three dimensions injection molded componentdesigns, it is generally desirable to start with a computerized solidmodeling package, such as Pro-Engineer™, CATIA™, I-DEAS™, Solid Works™,Solid Edge™, or other, which is used commonly in mechanical design anddrafting applications. The modeling package can be used to generatethree dimensional, photorealistic descriptions of the componentgeometry, called a solid model. At present, finite element analysiscodes based directly on solid models use nodes to define solid elementssuch as tetrahedra and hexahedra. In order to account for the physicsinvolved in injection molding, it is generally desirable to calculatefive quantities per node of the finite element model, namely pressure,three orthogonal components of velocity, and temperature. Given that asuitable model may contain hundreds of thousands of nodes, solution ofsuch complex numerical problems is difficult and requires substantialcomputer resources.

U.S. Pat. No. 5,835,379, issued to Nakano, and related European PatentApplication No. EP 0 698 467 A1, the disclosures of which areincorporated herein by reference in their entirety, suggest a method toreduce the number of variables to be determined in an injection moldingfinite element model during the filling phase in order to permitcalculation using lower computer resources than would be requiredotherwise. Nakano discusses the concept of flow conductance, K, toreduce the number of variables to two, namely pressure and flowconductance. As used herein, the terms flow conductance, fluidconductance, and fluidity are used interchangeably and are to beconsidered synonymous. The effect of varying material viscosity isincorporated into the calculation via the flow conductance variable.This necessarily involves the extrapolation of viscosity data and it iscontemplated that this leads to considerable error in the calculation ofviscosity. Also, during the filling and packing phases, the materialproximate the cavity walls begins to freeze and the thickness of thelayer increases with time until the part is ejected, which is thought tolead to additional error when applying the method of Nakano.

The viscosity, η, of a polymer melt is frequently measured as a functionof temperature and shear rate. When making a measurement, it is notpossible to measure the viscosity at low temperatures approaching thetemperature where the material solidifies, because at these lowtemperatures, the viscosity is relatively high and, at reasonable shearrates, we have discovered that thermal viscous dissipation issignificant. Accordingly, viscosity measurements are taken generally inthe range of temperatures at which the melt flows substantially readily.After measurement, a function is fitted to the measured data. Theavailable data fitting functions are generally reasonable at hightemperatures and hold over a wide range of shear; however, whenconsidering material at the wall that is below or near thesolidification temperature, it is necessary to extrapolate the viscositywell beyond the experimental range. This results in errors in theviscosity value, which in turn creates error in the flow conductance,because viscosity and flow conductance are related by the followingrelationship: $\begin{matrix}{{\nabla^{2}\kappa} = {- \frac{1}{\eta}}} & (1)\end{matrix}$which can also be represented in the form:∇·(η∇κ)=−1  (2)

Also, as noted by Nakano, the value of the flow conductance has a smallvalue at the cavity wall and increases with distance away from the wall.At the wall, the value of flow conductance is close to or equal to zero,because of the very high viscosity there. A common velocity boundarycondition in fluid flow is a zero velocity at the wall. This equates toa no-slip condition. When slip occurs, velocity at the wall is non-zero,and can be treated in several ways. Using the flow conductance approach,the flow conductance can be set to a small, non-zero value.Nevertheless, as a result, there are, in actuality, variations of flowconductance in the frozen layer which are orders of magnitude less thanin the melt. Such a wide variance in values leads to errors in thenumerical scheme of Nakano and others in calculating flow conductance.

During the packing phase, the values of pressure may also be obtained bysolving first for flow conductance. In packing, we have determined thatthe frozen layer becomes very significant, to the extent that eventuallythe majority of nodes are below the solidification temperature. Hence,the problems alluded to above are exaggerated in the packing phase.

SUMMARY OF THE INVENTION

Conventional three dimensional modeling techniques rely primarily on theunderlying laws of conservation of mass and conservation of momentum topredict fluid flow in an injection mold cavity. However, in order toachieve high accuracy and high predictability in a full, threedimensional simulation of the injection molding process, we havedetermined conservation of energy principles should be addressed aswell. In order to satisfy all three requirements simultaneously, withoutconsuming inordinate computational resources, creative methods have beendeveloped in accordance with the invention, based upon a variety of newmodeling assumptions and methodologies.

Accordingly, it is one object of the invention to employ conservation ofenergy principles in modeling fluid flow in a cavity. The fluid may be aviscous liquid, such as a polymer melt in an injection molding process,or a molten or semi-solid metal used in a casting process. Because thematerial also experiences a phase change from the liquid state to asolid state, considering the thermal effects underlying the conservationof energy principles both during initial flow and solidification havebeen found to be important.

For example, it has been recognized that during the filling and packingphases, there are three primary heat transfer mechanisms to consider:convection from the incoming melt, conduction out to the mold wall, andviscous dissipation, which is related to the thermal energy produced byshearing within the flowing polymer. Additionally, there may be othermechanisms, such as compressive heating effects, due to heat generatedby compression and cooling resulting from decompression. All threeprimary mechanisms contribute significantly to the energy balance duringfilling. Accordingly, accurate simulation of injection molding requirestypically non-isothermal analysis of the molten polymer as it flows intothe mold. During the packing phase, however, the flow of material isreduced significantly and the main heat transfer mechanisms areconduction to the mold wall and convection of melt from regions of highpressure to regions of lower pressure. Convection functions quitestrongly in the earliest stages of packing, when pressure is firstapplied. Note also that packing pressures can be lower or higher thanthe fill pressure, and can be profiled, that is they can be made tochange with time. Consequently material can flow into the mold and alsoout of the mold via the feed system, depending on the difference betweenthe pressure in the cavity and the applied pressure at the feed system.Analysis of the packing phase also typically requires that the melt beconsidered as a compressible fluid.

While conventional methods determine flow conductance at a plurality ofsmall elements that divide the region in which fluid flows, according toone embodiment of the present invention, flow conductance is calculatedat each of the nodes of each element. This technique allows a moreprecise value of the flow conductance to be defined, as the flowconductance may vary in any of a variety of desired fashions from nodeto node. The variation will depend on the type of finite element used todiscretize the problem. Use of a nodal value is more desirable than useof an elemental value, as other quantities such as temperature,pressure, and velocity are more precisely defined as point values and,again, may be interpolated from node to node in any fashion. Thisinterpolation capability also enables the location of the interfacebetween the solid and molten polymer to be accurately determined, whichhas been found to be particularly relevant in simulation predictiveaccuracy.

Accordingly, it is an object of the invention to use a one dimensionalanalytic function to describe the local temperature distribution at anode. It is another object of the invention to define the variation of aone dimensional analytic function, with time, to account for heatconvection. It is yet another object of the invention to define thevariation of a one dimensional analytic function to account for viscousheat generation. It is still another object of the invention to definean explicit temperature convection scheme using the one dimensionalanalytic function. And it is a further object of the invention toprovide anisotropic finite element mesh refinement, including a distanceto the cavity wall computation algorithm.

Moreover, according to one embodiment, the present inventionincorporates the concept of a frozen layer to remedy certain of theproblems with respect to conventional simulation techniques mentionedabove. A criterion is used to define the interface between the moltenpolymer and the solidified polymer at the mold cavity wall. Thiscriterion may be a temperature, a minimum value of velocity, acombination of the two, or some other physical quantity. In the frozenregion, there is no significant movement of polymer and so it is notnecessary to calculate the flow conductance there. Nodes that aredetermined to be in the frozen region can be removed from the solutionprocess for velocity, pressure, and fluid conductance. Only temperatureneed be calculated in the frozen layer. Note, however, that propertiesother than temperature can also be computed in the frozen layer. Forexample, the material morphology will change with change in state.Crystallinity of semi-crystalline polymers changes with temperature, andcrystals can only be said to form once the melt has cooled sufficientlyfrom the melt state. The state of stress also depends on temperature anddiffers between liquid and solid phases. Polymers are viscoelastic, sotheir stress state, even in the solid, changes with time andtemperature.

In effect, by only calculating temperature, this reduces the number ofvariables to be determined to one in the frozen layer and therebydramatically increases the speed of computation, while reducing thememory required. By removing the frozen nodes, the need forextrapolation of the viscosity is removed and so the errors associatedtherewith are removed as well. Further, by removing nodes in the frozenlayer from the analysis, very small values of flow conductance areremoved and the equations to be solved are better conditioned than inconventional methods. Still further, by not calculating fluidconductance in these regions, it has been found that not only is theproblem size reduced considerably, but the stability of numericalmethods is improved as well.

Accordingly, it is another object of the invention to determine theposition of the solid/liquid interface in the mold cavity and withinelements defined in the mesh. It is still another object of theinvention to develop a formulation of elements with a linear variationof material properties throughout. It is yet another object of theinvention to determine an effective viscosity function in elementscontaining the solid/liquid interface. It is a further object of theinvention to eliminate frozen nodes and elements from the solutiondomain. And it is yet a further object of the invention to compute theeffective pressure in regions that have solidified.

One embodiment of the invention involves a method for modeling injectionof a fluid into a mold defining a three dimensional cavity, the methodincluding the steps of providing a three dimensional solid computermodel defining the cavity, discretizing a solution domain based on thesolid model, specifying boundary conditions, solving for filling phaseprocess variables using conservation of mass, momentum, and energy toprovide respective filling solutions for at least the portion of thesolution domain, then solving for packing phase process variables in asimilar manner to provide respective packing phase solutions for atleast the portion of the solution domain, and determining whether thefilling and packing phase solutions are acceptable. The filling andpacking phase process variables can include density, fluidity, moldcavity fill time, mold cavity packing time, pressure, shear rate, shearstress, temperature, velocity, viscosity, and volumetric shrinkage. Inthe event the filling and/or packing phase solutions are unacceptablethe boundary conditions and or discretized solution domain can bemodified and the analysis repeated until an acceptable result isachieved. To facilitate a user in determining whether the results areacceptable, a variety of filling and packing phase solutions can bedisplayed graphically, such as fill time, packing time, density,pressure, shear rate, shear stress, temperature, velocity, viscosity,and volumetric shrinkage.

Another embodiment of the invention involves a method for modelinginjection of a fluid into a mold defining a three dimensional cavityincluding the steps of providing a three dimensional solid computermodel defining the cavity, discretizing a solution domain based on thesolid model, specifying boundary conditions, solving for filling phaseprocess variables using conservation of mass, momentum, and energy toprovide filling phase solutions for at least some of the portion of thesolution domain, and determining whether the solutions are acceptablefor injection of the fluid during filling of the mold cavity.

The discretizing step may include the substep of generating a finiteelement mesh based on the solid model by subdividing the model into aplurality of connected elements defined by a plurality of nodes. Themesh generating substep may include generating an anisotropic mesh inthick and thin zones of the model, such that mesh refinement providesincreased resolution in a thickness direction without increasingsubstantially mesh refinement in a longitudinal direction.

The boundary conditions may include such parameters as fluidcomposition, fluid injection location, fluid injection temperature,fluid injection pressure, fluid injection volumetric flow rate, moldtemperature, cavity dimensions, cavity configuration, and mold partingplane, and variations thereof.

Solving for the filling phase process variables may include the substepsof solving for fluidity, pressure, velocity, and viscosity for at leastsome of the solution domain, where viscosity is based on temperature.Temperature, in turn, may be based on convective, conductive, and/orviscous dissipation heat transfer contributions. Velocity and/orviscosity may be calculated iteratively, until pressure converges.

This method may also include the substep of determining free surfaceevolution of the fluid in the cavity based on velocity in a given timeincrement, wherein the free surface evolution is determined iteratively,until the cavity is filled.

Once the simulation has modeled that the cavity is filled, the methodmay then solve for packing phase process variables using conservation ofmass, momentum, and energy for at least a portion of the solution domainbased on respective states of the process variables at termination offilling, to provide respective packing phase solutions for at least someof the solution domain. A determination can then be made whether thepacking phase solutions are acceptable for injection of the fluid duringpacking of the mold cavity.

Again, solving for packing phase process variables can include solvingfor fluidity, pressure, velocity, and viscosity for at least some of thesolution domain, where viscosity is based on temperature. Temperature,in turn, may be based on convective, conductive, and/or viscousdissipation heat transfer contributions. Velocity and/or viscosity maybe calculated iteratively, until pressure converges. The mass propertiesof a component produced in accordance with the boundary conditions canalso be determined. The mass properties may include component density,volumetric shrinkage, component mass, and component volume and may alsobe calculated iteratively, along with velocity and viscosity, until apredetermined pressure profile is completed.

BRIEF DESCRIPTION OF THE DRAWINGS

The novel features believed characteristic of the invention are setforth and differentiated in the appended claims. The invention, inaccordance with preferred and exemplary embodiments, together withfurther advantages thereof, is more particularly described in thefollowing detailed description taken in conjunction with theaccompanying drawings in which:

FIG. 1 is a schematic representation of computer hardware apparatussuitable for use with the disclosed methods for modeling injection offluid in a mold cavity in accordance with one embodiment of the presentinvention;

FIG. 2 is a schematic representation of a top level system flowchartsummarizing certain process steps in accordance with one embodiment ofthe present invention;

FIGS. 3A and 3B are schematic representations of an interior point inthe mold cavity and the related temperature distribution along a linejoining the interior point and the mold wall;

FIG. 4 is a schematic representation of the association of a onedimensional analytic function for temperature distribution at eachpoint, i, within a mold cavity;

FIG. 5 is a schematic representation in flowchart form of the stepsinvolved in an overall temperature solution over a global time incrementin accordance with one embodiment of the present invention;

FIG. 6 is a schematic representation of velocity at a node and therelationship thereof with surrounding elements;

FIG. 7 is a schematic representation of velocity on a face of anelement;

FIG. 8 is a schematic representation of velocity at a node and therelationship thereof with a connected, upstream tetrahedral element;

FIG. 9 is a schematic representation of velocity at a face of an elementand the relationship thereof with a connected, upstream tetrahedralelement;

FIGS. 10A, 10B are schematic representations of an interior point in themold cavity and interpolation of upstream temperature using a onedimensional analytic function at upstream element nodes;

FIG. 11 is a schematic representation of tetrahedral element facesplitting combinations;

FIG. 12 is a schematic representation of a template for splitting atetrahedral element on one edge;

FIG. 13A is a schematic representation of a template for splitting atetrahedral element on two adjacent edges;

FIG. 13B is a schematic representation of a template for splitting afour-sided base pyramid;

FIG. 13C is a schematic representation of a template for splitting atetrahedral element on two opposite edges;

FIG. 14A is a schematic representation of a template for splitting atetrahedral element on three edges with a common shared face;

FIG. 14B is a schematic representation for splitting a tetrahedralelement on three edges with a common node;

FIG. 14C is a schematic representation of a template for splitting atriangular prism;

FIG. 14D is a schematic representation of a template for splitting atetrahedral element on three edges in series;

FIG. 15A is a schematic representation of a template for splitting atetrahedral element on four opposite edges;

FIG. 15B is a schematic representation of a template for splitting atetrahedral element on four adjacent edges;

FIG. 16 is a schematic representation of a template for splitting atetrahedral element on five edges;

FIG. 17 is a schematic representation of a template for splitting atetrahedral element on six edges;

FIGS. 18A and 18B are schematic representations of node layer numberingand connectivity trees used in determination of distance to the cavitywall to limit the possible search space;

FIG. 19 is a schematic representation of a method for determining thesolid/liquid interface position within an element using solidificationtemperature and a one dimensional analytic function for temperaturedistribution at each liquid node;

FIG. 20 is a schematic representation of the variation of viscositywithin an element which contains a solid/liquid interface;

FIG. 21 is a schematic representation of fluid flow in the vicinity of asolid/liquid interface;

FIG. 22 is a schematic representation of a method for determining a corenode associated with each wall node and each internal node;

FIGS. 23A through 23C are schematic representations of the relationshipbetween wall, internal, and core nodes, and a method for determiningpressure at frozen nodes;

FIG. 24 is a schematic representation of a filling phase flowchartsummarizing certain process steps in accordance with one embodiment ofthe present invention; and

FIG. 25 is a schematic representation of a packing phase flowchartsummarizing certain process steps in accordance with one embodiment ofthe present invention.

DETAILED DESCRIPTION OF THE INVENTION

Modeling of the injection molding process in a three dimensionalsimulation can be described by the conservation equations for mass,momentum, and energy, respectively, as follows: $\begin{matrix}{{\frac{\partial\rho}{\partial t} + {\nabla{\cdot \left( {\rho\quad\overset{\rightarrow}{v}} \right)}}} = 0} & (3)\end{matrix}$where ρ is density, v is velocity, t is time, ∇ stands for the gradientwith respect to a position vector, and → denotes a vector quantity.$\begin{matrix}{{{\rho\left( {\frac{\partial\overset{\rightarrow}{v}}{\partial t} + {\overset{\rightarrow}{v} \cdot {\nabla\overset{\rightarrow}{v}}}} \right)} + {\nabla{\cdot P}}} = {\rho\quad\overset{\rightarrow}{g}}} & (4)\end{matrix}$where P is the momentum flux and g is the gravity force. $\begin{matrix}{{{\rho\left( {\frac{\partial\hat{U}}{\partial t} + {\overset{->}{v} \cdot {\nabla\hat{U}}}} \right)} + {P:{D + {\nabla{\cdot \overset{->}{q}}}}}} = {\rho\quad Q}} & (5)\end{matrix}$where U is specific internal energy, D is the deformation rate tensor, qis heat flux, and Q is the specific heat rate which signifies other heatsources or sinks.

Equation 4 is commonly expressed in terms of temperature and pressure,and for the case in which Q=0, i.e. no additional heat sources or sinks,typical for thermoplastics, takes the form, $\begin{matrix}{{{\rho\quad{C_{p}\left( {\frac{\partial T}{\partial t} + {\overset{->}{v} \cdot {\nabla T}}} \right)}} + {\frac{T}{\rho}\frac{\partial\rho}{\partial T}\left( {\frac{\partial p}{\partial t} + {\overset{->}{v} \cdot {\nabla p}}} \right)}} = {{\eta\quad{\overset{.}{\gamma}}^{2}} + {\nabla{\cdot \left( {\kappa{\nabla T}} \right)}}}} & (6)\end{matrix}$where C_(p) is specific heat, T is temperature, p is pressure, η isviscosity, is shear rate, and κ is thermal conductivity. This κ shouldnot be confused with flow conductance.

In general, these equations are analytically intractable, so they aresolved by a discrete numerical procedure, such as a boundary elementmethod, a finite element method, a finite difference method, a finitevolume method, or a meshless method. According to the finite elementmethod, solution requires subdivision of the solution domain into a setof smaller sub-domains, and transformation of the equations into a setof discrete equations describing the ensemble of the smallersub-domains. Application of the finite element method to a threedimensional domain results in three dimensional sub-domains calledelements. According to one method, four node, linear tetrahedralelements may be used, wherein the elements are defined by nodes at thevertices of each tetrahedral element, and the term linear refers to thepolynomial order of the element interpolation or shape functions.

Most common finite element practice results in a discretization of thedomain, such that the field variables of interest are computed at thenode points, and the material properties required for analysis areconstant within an element. To improve accuracy, especially in highgradient regions, it has been determined to be desirable to formulatethe equations such that both the field variables (such as pressure,temperature, and velocity) and the material properties (such asviscosity) are defined at the nodes and are, unless otherwise stated,interpolated within an element according to the element shape function.Discretization in this way has been determined to provide not onlyadditional computational accuracy, but also facilitates use of higherorder advanced interpolations, to be described more fully hereinbelow.

Ultimately, the accuracy of simulation of a physical process dependsstrongly on how well the model used in the simulation can describeimportant phenomena that are characteristic of the process. Forinjection molding, there are a number of specific phenomena that reflectthe complex nature of the process, and which need to be modeled well foraccurate predictive simulation results. Among these phenomena are: verylarge and strongly varying temperature gradients near the solid/liquidinterface; variation in position of the solid/liquid interface boundarywith time and location in the mold cavity; large variations in polymermaterial properties at all places and at all times during molding; andincreasing volume of solid material with time as the injected liquidmaterial solidifies during molding.

It is possible to capture these specific phenomena in a number of ways.One method is to refine the mesh. This equates to reducing the elementsize sufficiently to ensure that linear interpolations withinneighboring elements can adequately capture rapid spatial changes in thefield variables. However, a downside to this approach is that the numberof elements and nodes increases enormously, often making thecomputational cost prohibitive. This applies in particular to injectionmolding simulation where part dimensions such as thickness and lengthcan differ by two orders of magnitude or more and where the largestspatial gradients of the field variables are found across the smallestspatial dimension, namely the thickness.

An alternative is to use elements with higher order shape functions,such as polynomial functions. Higher order shape functions enablenon-linear interpolations of the field variables within an element, withthe prospect that fewer elements are required. However, the number ofnodes required to define the higher order shape functions within theseelements increases and consequently the system of equations is stillprohibitively large.

Accordingly, it has been determined that certain modeling functions andassumptions can be employed alone, or in combination, to significantlyenhance simulation accuracy, without increasing the number of nodessubstantially, with a concomitant detrimental impact on computation timeor resources.

Insofar as simulation of the injection molding process entailssignificant computational effort, depicted in FIG. 1 is a computerhardware apparatus 10 suitable for use with the disclosed methods formodeling injection of fluid in a mold cavity in accordance with thisinvention. The apparatus 10 may be a portable computer, mini-computer,or other suitable computer having the necessary computational speed andcapacity to support the functionality discussed more fully hereinbelow.The computer 10 typically includes one or more central processing units12 for executing the instructions contained in the software code whichembraces the simulation model. Storage 14, such as random access memoryand read only memory, is provided for retaining the model code, eithertemporarily or permanently, as well as other operating software requiredby the computer 10. Permanent, non-volatile read/write memory such ashard disks are typically used to store the code, both during its use andidle time, and to store data generated by the software. The computer 10also includes one or more inputs 16, such as a keyboard and disk readerfor receiving input such as data and instructions from a user, and oneor more outputs 18, such as a monitor or printer for providingsimulation results in graphical and other formats. Additionally,communication buses and I/O ports are provide to link all of thecomponents together and permit communication with other computers andcomputer networks, as desired.

FIG. 2 is a schematic representation of a highly simplified, top levelsystem flowchart summarizing certain process steps in accordance withone embodiment of the present invention. As a first step 20, a threedimensional solid computer aided design model of the mold cavity isgenerated or provided, as discussed above. The model solution domain isthen defined and discretized by any of a variety of methods, such as byfinite element analysis in which a finite element model is produced bygenerating a finite element mesh based on the solid model in step 30.The mesh consists of a plurality of contiguous solid elements defined byshared nodes. With the resultant finite element model or otherdiscretized solution domain defined, a user specifies boundaryconditions in step 40 for the analysis. These boundary conditions, byway of example, may include geometric constraints and numerical valuesrelated to fluid composition, fluid injection location, fluid injectiontemperature, fluid injection pressure, fluid injection volumetric flowrate, mold temperature, cavity dimensions, cavity configuration, moldparting plane, and other initial design data. Moreover, these conditionsmay not be constant values, but rather time dependent or distributedprofiles. For example, injection pressure may vary with time or the moldtemperature may be different in different zones, as a result of a priorcooling analysis or active cooling.

Once the boundary conditions have been entered, the computer 10 executesthe instructions in accordance with the simulation model to firstcalculate or solve relevant filling phase process variables for thenodes in step 50. As will be discussed in great detail hereinbelow withrespect to the filling phase flowchart depicted in FIG. 24, suchvariables can include fluidity, mold cavity fill time, pressure, shearrate, stress, velocity, viscosity, and temperature. Also what can becalculated is not limited to these variables; however these are basicvariables that can be used to solve other variables included incalculations of such things as crystallization kinetics and fiberorientation distributions. Further, filling can also be solved as acompressible fluid, in which case mass terms included in the packingphase calculations (e.g. density, mass, and volumetric shrinkage) canalso be calculated in the filling phase. According to one embodiment,the simulation can be based on the assumption that the fluid isincompressible in the filling phase and compressible in the packingphase. According to another embodiment, it can be assumed that the fluidis compressible in both the filling and packing phases. Also, as will bediscussed in further detail below with respect to FIGS. 24 and 25, it isnot mandatory to solve for fluidity before pressure, velocity, andviscosity, nor is it necessary to solve for fluidity at all. This isgenerally only required when a method such as that discussed in Nakanois used to solve the conservation of mass and momentum equations.Conservation of mass and momentum can also be solved using formulationssuch as Stokes and Navier-Stokes. Accordingly, the principles presentedherein have application to all of these descriptions of the physicalprocess such as Nakano, Stokes, and Navier-Stokes.

Once the simulation reaches the stage in the analysis where it isdetermined that the mold cavity has been filled, the computer 10executes the instructions in accordance with the simulation model tonext calculate or solve relevant packing phase process variables for thenodes in step 60. As will be discussed in great detail hereinbelow withrespect to the packing phase flowchart depicted in FIG. 25, suchvariables can include the mass properties of the component produced inaccordance with the simulation model such as density and volumetricshrinkage, in addition to fluidity, packing time, pressure, shear rate,stress, velocity, viscosity, and temperature.

Upon completion of the analysis, the analytical results may be output instep 70 in a variety of manners. For example, the relevant variables maybe displayed in a graphics format, overlaying the solid model for visualreview by the user, or may be output electronically for furtherprocessing or analysis. If the results of the filling phase and thepacking phase are deemed to be acceptable in step 80, the simulationterminates in step 100 and the user can proceed to release the design tomanufacturing. Because the specified boundary conditions includedinformation related to the configuration of the injection mold and theprocess parameters, the design can be released for machining of theinjection mold and the injection molding process operation sheetsgenerated directly.

If, however, the user determines that the results of the simulation instep 70 are unacceptable or less than optimal, the user has the optionin step 90 of modifying one or more of the boundary conditions and/ordiscretization of the model solution domain and thereafter repeatingsimulation steps 50 through 70 iteratively, until such time as the useris satisfied with the results. Examples of unacceptable results includeanalytical instability of the model or process failures such as shortshots, wherein the mold cavity is incompletely filled, or generation ofexcessive temperatures, velocities, or pressures during filling whichcould degrade component polymer material properties or introduceexcessive residual stresses in the components which would adverselyeffect production yields and could lead to premature component failure.By providing this highly accurate analytical simulation capability earlyin the design process, significant costs and delays downstream duringinitial production runs can be avoided.

Before looking at the process steps of the filling phase and packingphase flowcharts of FIGS. 24 and 25, it may be advantageous to developthe underlying theory and basis for certain of the approaches takentherein.

As discussed hereinabove, certain limitations of conventional threedimensional injection molding modeling techniques relate to the largetemperature gradients, movement of the solid/liquid interface boundary,large material property variations during solidification, andsolidification of melt material during molding. Loosely grouped, aboutfive innovative methodologies related to advantageous application ofconservation of energy principles have been identified which offersubstantial improvements in predictive accuracy and/or computation time.The first is the use of a one dimensional analytic function to describethe local temperature distribution at a node. The next is the definitionof the variation of a one dimensional analytic function, with time, toaccount for heat convection. The third relates to the description of thevariation of a one dimensional analytic function to account for viscousheat generation. The next entails the description of an explicittemperature convection scheme using the one dimensional analyticfunction. And the fifth addresses finite element mesh refinement,including a distance to the cavity wall computation algorithm. Whilerecognition and application of the conservation of energy principles tothe modeling of injection molding simulation are being categorized herein five general areas, these are not to be construed as limiting orotherwise, and are being presented in this manner solely for purposes ofclarity.

Looking first to use of a one dimensional analytic function to describethe local temperature distribution at a node, at any point in timeduring mold filling and packing, a temperature distribution existsbetween any given interior point within the material in the cavity andthat part of the mold wall, in the immediate vicinity of the interiorpoint. In general, the temperature distribution can be taken along aline that is the shortest distance, d, between the interior point andthe mold wall. See, for example, FIG. 3A. FIG. 3B is a schematicrepresentation of the temperature profile in the fluid along a linejoining an interior point in the mold cavity and the mold wall, as afunction of distance, d, from the mold wall.

A variety of analytic functions or discrete functions can be used todescribe the temperature distribution. According to one embodiment, anerror function formulation can be used, because of its physicalsignificance in the solution of conduction heat transfer in asemi-infinite solid. According to the invention, any and every interiorpoint is treated as if it were part of the domain of a semi-infinitesolid. Each point is treated independently of every other point. Theequation for heat conduction in a semi-infinite solid is used, withmodification, to provide a better estimate of the local temperaturedistribution at a point in the cavity at all locations and times duringthe injection molding process, and respecting all significant heattransfer mechanisms.

The temperature distribution in a semi-infinite solid varies with timeand distance from the finite boundary. In this application, the finiteboundary may be defined as the mold cavity wall. When the temperature atthe finite boundary of a semi-infinite solid, initially at a uniformtemperature, T_(∞), is instantly changed to another temperature, T_(w),the distribution of temperature in the solid varies with distance fromthe wall and time according to the following function:T=T _(w)+(T _(∞) −T _(w))erf(d/ 2{square root}{square root over(αt)})  (7)where:

-   erf=error function;-   d=distance from wall;-   t=time;-   T_(∞)=temperature at infinite distance from the wall;-   T_(w)=wall temperature; and-   α=thermal diffusivity of solid.

See, for example, Transport Phenomena, by R. B. Bird, W. E. Stewart, andE. N. Lightfoot, published by John Wiley & Sons, 1960.

This function has been found to have several useful properties in theapplication to injection molding simulation. For example, it can capturespatial and temporal variation in temperature at both short and longdistances and times. It is a simple form and easily computed. There arefew adjustable parameters. It has physical significance, since heatconduction is a major heat transfer mechanism, and at short contacttimes, the equation is a close approximation to the solution for afinite slab, that is for a solid bounded by walls on both sides, similarto the situation in a mold cavity.

This function may be used advantageously to improve the localtemperature accuracy, because over short time periods, the distributionof temperature is highly non-linear over a relatively large distanceduring injection molding. Since it is generally desirable to keep thenumber of nodes and elements relatively low in order to reduce thecomputational cost, the element dimensions are relatively large comparedto the local cavity geometry being modeled. Also, the elements havelinear shape functions. Consequently they cannot adequately describe thevariation of temperature under circumstances where the temperaturedistribution is highly non-linear. This error function formulation candescribe such a non-linear distribution of temperature.

In its application, this function may be defined for each and everyinterior node in the three dimensional domain. The equation parametersfor each node are defined separately for each node. See, for example,FIG. 4, which shows schematically the association of a one dimensionalanalytic function for temperature distribution at each point, i, withina mold cavity. Three points are shown.

The next innovative methodology relates to the definition of thevariation of a one dimensional analytic function, with time, to accountfor heat convection. The use of the error function equation in the formpresented hereinabove in certain circumstances may not be sufficient tomodel adequately the unique variation in temperature of a fluid duringinjection molding. In a strict sense, the error function equation ismeant to be applied to a fluid that is semi-infinite in extent, andstationary. However, in a real injection molding application, the fluidis finite in extent, being bounded by the cavity walls, and the fluid isin motion for much of the process. The fluid motion, in particular,introduces another significant heat transfer mechanism, namelyconvection, the contribution of which should be included in the overallconservation of energy analysis to improve simulation predictiveaccuracy.

Consider the physical process of injection molding. Material contactenables heat conduction to proceed, which, in the absence of other heattransfer mechanisms, is adequately described by the error functionequation at short times. However, convective heat transfer issuperimposed on this mechanism. Convective heat transfer has the effectof retarding heat loss by conduction, because injection moldingcontinuously brings hot material from upstream to any point underconsideration during the filling phase.

Notice that the error function equation includes a term for time.Accordingly, it has been determined that one way to address the thermalinput is to consider that the convective transfer of heat has the effectof retarding the passage of time for conduction at any point underconsideration. So whilst real global clock time proceeds at a constantrate, describing the passage of injection time from filling throughpacking, each node may be considered to carry its own node thermalclock, which varies locally and typically may not proceed at the samerate as the global clock.

For example, if the fluid at an interior point in the cavity isstationary, then for as long as this is true, the node thermal clockmarks out time at the same rate as the global clock. If, however, thefluid at an interior node is moving very quickly so that convection heattransfer dominates, then any local heat loss by conduction is whollycounteracted by heat input due to convection. In this case, the nodethermal clock is stationary, not advancing at all, even while the globalclock has moved on.

In the case where the flow is intermediate, so that neither conductionnor convection dominate, then the node thermal clock for an arbitrarynode, i, may be incremented according the following equation: (8)t _(i) ^(k) =t _(i) ^(k−1)+exp(−F_(pe) Pe _(i))Δt  (8)where:

-   t_(i) ^(k)=node thermal clock time used in the error function    equation for node i at time step k;-   t_(i) ⁻¹=node thermal clock time used in the error function equation    for node i at a previous time step k−1;-   Pe_(i)=Peclet number for node i;-   F_(pe)=tuning parameter; and-   Δt=global time clock increment.    and further where:    Pe _(i) =d _(i) ·v _(i)/α_(i)  (9)    where:-   d_(i)=distance from wall of node i;-   v_(i)=local velocity of node i; and-   α_(i)=thermal diffusivity of node i.

Whilst any of a number of equations could be used to describe the nodethermal clock time advance, this empirical equation was selected becauseit has the desired properties. In particular:As σ→0, exp(−F_(pe)Pe)Δt→Δt;As σ→∞, exp(−F_(pe)Pe)Δt→0;As α→0, exp(−F_(pe)Pe)Δt→0; andAs α→∞, exp(−F_(p)ePe)Δt→Δt.

The tuning parameter constant F_(pe) is found by tuning simulationresults to experimental data. The data can be explicit temperaturedistribution data from injection molding experiments or indirect data,such as injection pressure measurements. It has been determined that theconstant F_(pe) may generally be set in the range of between about zeroand one.

Yet another innovative methodology addresses the variation of a onedimensional analytic function to account for viscous heat generation.Since movement of the fluid requires mechanical work, thermodynamicprinciples require that heat is also added to the system through thiswork. In the case of fluid motion, this work takes the form of viscousdissipation or frictional heat. Frictional heat varies with timethroughout the interior of the mold cavity during injection molding. Inany time period, this thermal energy contribution is given by theproduct of the local shear stress and the local shear rate over thattime period. It has been determined that this thermal contribution tendsto be most significant in the vicinity of solid/fluid boundaries, whereboth viscosity is high, due to the relatively low temperature proximatethe frozen material, and stress is high, due to motion of this viscousfluid.

Frictional heat can be considered as a local, time variant heat sourceat each point in the cavity. This local addition of heat cannot beaccounted for by the above-referenced error function equation in itsstandard form. Up to this point, there is no mechanism to increase thefluid temperature in the functional description above the maximum valuedefined by the value of T_(∞); whereas, in injection trials, it has beendetermined that the temperature of the fluid proximate the liquid/solidboundary can exceed not only T_(∞), but also the temperature of thefluid at the point of injection into the mold cavity.

Accordingly, to accurately model the temperature of the fluid in thecavity in a simulation, frictional heat should be considered in theanalysis. When this is done, it may occur that the local temperature,even after convection and conduction, has increased above T_(∞). Thisshift in temperature can be accommodated by adjusting the far fieldtemperature, T_(∞), at each time step. The adjustment is done byrecalculating T_(∞) at the end of each time step in the simulation, forexample, by substituting the computed nodal temperature after solutionof the energy equation, into the error function equation and deriving anew value of T_(∞) at each node, i, in the cavity as follows:$\begin{matrix}{T_{\infty_{i}} = {T_{w_{i}} + \frac{\left( {T_{i} - T_{w_{i}}} \right)}{{erf}\quad\left( {{d_{i}/2}\sqrt{\alpha_{i}t_{i}}} \right)}}} & (10)\end{matrix}$

For any interior node, therefore, there is a mechanism to describe thelocal one dimensional temperature distribution and its variation inspace and time. This permits the use of a better approximation than asimple linear interpolation, or even a higher order polynomialinterpolation afforded by the element shape function. The new functioncan be used to dramatically improve the accuracy of convective heattransfer calculations and to compute the position of the solid/liquidinterface at any time during the injection molding process.

FIG. 5 is a schematic representation in flowchart form of the stepsinvolved in the overall temperature solution over a global timeincrement, Δt, according to one embodiment of the invention, the detailsof which will be discussed more fully hereinbelow. In a first step 15,the contribution of convective heat transfer over the time increment iscalculated. Then in step 25, the contribution of viscous dissipationheat transfer is calculated. Thereafter, in step 35, the total energyequation is solved, taking into account heat conduction over the timeincrement. The temperature T_(∞) is then updated, for each node, and theresults are available for use in solving the momentum equation,determining evolution of the flow front during the filling phase, andquantifying material property variation, as will be discussed in greaterdetail hereinbelow.

A further area which has been addressed by the current inventionconcerns the description of an explicit temperature convection schemeusing the one dimensional analytic function. Solution of the energyconservation equation for injection molding requires, in part, thecomputation of the contribution of convection to the overall heattransfer process. While a number of mechanisms have been contemplated todeal with convective heat transfer, according to one embodiment of theinvention, an explicit scheme is used, in which the contribution totemperature change at each node is computed individually.

Since an explicit scheme entails finding the upstream point which willflow into a target node in a given time increment and interpolation ofthe temperature at that upstream point, the previously described onedimensional analytic function can be used to provide a far more accurateinterpolation of temperature at the upstream point than otherwise. ALagrangian interpolation is employed to determine the upstream terms,otherwise known as the advection terms. The same mechanism can be usedto solve both the advection of temperature (convection), required forthe energy equation solution, and the advection of matter(concentration) required for the solution of the free surface or flowfront position within the mold cavity during the filling phase.

The Peclet numbers for the energy and concentration function equationsare very large and infinity, respectively. Conventional finite elementmethods to simulate this feature are either unstable or tend toincorporate an unacceptable degree of diffusion that destroys accuracy.This is because conventional least-squares computational methods areunnatural ways to deal with the problem. It has been determined that thenatural way to deal with hyperbolic differential equations is to followthe function along its characteristics. In the case of fluid flow in theinjection molding context, this translates to tracking the path of anarbitrary “particle” in the continuum.

One object of the algorithm is, given any point in a three dimensionaldomain (in particular, a node), to use the velocity field to track theorigin of the particle that now occupies the node, over the last timestep. In other words, the concept is to look upstream to see where theparticle, which now occupies the node, was at one time increment ago.Once the origin is determined, an interpolation can be performed todetermine the value of the concentration function or temperature there,at the upstream point. Then, given this value, the value at the nodeunder consideration can be updated accurately in the present timeincrement.

With reference now to FIG. 6, starting with a node, the velocity vectoris as shown at node A. Without loss of generality, the elements in thefigures are depicted as triangles, rather than tetrahedra, so as tofacilitate depiction; however, as stated hereinabove, the elements mayalso comprise hexahedra, or other suitable elements shapes, depending onthe mesh method utilized and the configuration of the model.

A cursory review of the figure suggests that the particle currentlyoccupying node A comes from the triangle that has the opposite face L.However, an algorithm recognizes this by testing to see whichtriangle/tetrahedron satisfies all required criteria. The test should becomputationally fast. This part of the calculations, of determiningwhich of the surrounding triangles/tetrahedra is the correct one, isperformed in subroutine upwelm1.F, discussed hereinbelow.

Once the correct triangle/tetrahedron is recorded, it may be that theparticle comes from beyond the face L. This may be confirmed bycomparing the current time step size with the time required for aparticle to traverse the triangle/tetrahedron starting from face L. Ifthe search is required to extend beyond L into a contiguoustriangle/tetrahedron, subroutine contig.F, also discussed hereinbelow,determines the correct neighboring triangle/tetrahedron and rearrangesthe local shape function values accordingly. In a triangle example,suppose the shape function array at L is [0.2,0.8,0.0]. Then if theorder of the node numbers on L are reversed in the contiguous element,the shape function array becomes [0.8,0.2,0.0], it being assumed thatthe third node and A both are the last nodes in the elements, ergo thevalue of 0.0 at face L.

The next issue is represented schematically in FIG. 7. The velocity iscomputed at the relevant position on L, and is used to determine fromwhich face a particle with that velocity could originate. This iscomputed in subroutine upwelm2.F, discussed below. One will note thatthe main difference between upwelm1.F and upwelm2.F is that in theformer the search is from a node for the upstream element of theconnected elements; whereas, in the latter, the search is from a facefor the upstream face within the same element. Also note, that if theparticle comes from beyond the opposite face, then after rearranging thenumbering in contig.F, all further searches can be accomplished inupwelm2.F.

Both upwelm1.F and upwelm2.F subroutines start from the fundamentalrelation between the global and local coordinates:x−x ₄ =J·L  (11)The letters are displayed in bold to denote a vector (x, L) or a matrix(J), where: x is any global position vector (coordinates); L is thelocal coordinate vector; J is the Jacobian of transformation between thetwo coordinate systems; and x₄ is the location of the fourth node of theelement in the global coordinate system.

The upwelm1.F subroutine may be described as follows, with reference toFIG. 8. Starting at a node A, it is desirable to find the first upwindintercept with an element facet, or the upwind location within theelement if the time-step is not sufficient to reach the facet. Thevelocity at A is known, as depicted in FIG. 8.

In a first case, suppose node A is the fourth node in the element nodeorder of the tetrahedron shown in the figure. This yields:L=J ⁻¹*(x−x ₄)  (12)

The next step is to find where the particle trajectory intercepts the(infinite) plane defined by the facet opposite the node A. Note thatthere must be an intercept somewhere, unless the velocity is parallel tothe plane. A position x is on the trajectory of the particle ifx=x ₄ −kv  (13)Here, k is the time upstream from x₄. Substituting equation (13) inequation (12)12), the local coordinates of this position x are:L=J ⁻¹·(−kv)  (14)The trajectory intercepts the plane where L₄=0, therefore:$\begin{matrix}{{\sum\limits_{i = 1}^{3}L_{i}} = 1} & (15)\end{matrix}$Expanding this yields: $\begin{matrix}{{{- k}\quad{\sum\limits_{i = 1}^{3}{\sum\limits_{j = 1}^{3}{J_{ij}^{- 1}v_{j}}}}} = 1} & (16)\end{matrix}$Making k the subject of the formula: $\begin{matrix}{k = \frac{- 1}{\sum\limits_{i = 1}^{3}{\sum\limits_{j = 1}^{3}{J_{ij}^{- 1}v_{j}}}}} & (17)\end{matrix}$

For this facet to be a candidate, k must be both positive (i.e. a facetupstream of node A) and finite. An infinite value implies that thevelocity is parallel to the opposite face L. If k is both positive andfinite, the next step is to check if the intercept with the facet planeis within the boundaries of the facet. This is true if and only if allthree components of L associated with the facet nodes are non-negative.The components can then be computed by substituting k from equation (17)into equation (14).

In a second case, suppose A is not the fourth node, but some node m,where m is any of 1, 2, or 3. Then the intercept of the trajectory withthe plane of the facet opposite A occurs where L_(m)=0. Suppose this isat some position x so that:x=x _(m) −kv  (18)Inserting equation (18) in equation (12), and considering the mth node(the notation uses m to describe both local and global node number),yields: $\begin{matrix}{L_{m} = {{\sum\limits_{j = 1}^{3}{J_{mj}^{- 1}\left( {x_{m}^{j} - {k\quad v_{j}} - x_{4}^{j}} \right)}} = 0}} & (19)\end{matrix}$

The superscript on x indicates the coordinate component (for example,x²≡y, and x³≡z).Making k the subject of the formula yields: $\begin{matrix}{k = \frac{\sum\limits_{j = 1}^{3}{J_{mj}^{- 1}\left( {x_{m}^{j} - x_{4}^{j}} \right)}}{\sum\limits_{j = 1}^{3}{J_{mj}^{- 1}v_{j}}}} & (20)\end{matrix}$

For this facet to be a candidate, k must be both positive (i.e. a facetupstream of node A) and finite. An infinite value implies that thevelocity is parallel to the opposite face L. If k is both positive andfinite, the next step is to check if the intercept with the facet planeis within the boundaries of the facet. This is true if and only if allthree components of L associated with the facet nodes are non-negative.The components can then be computed by substituting k from equation (20)into equation (14).

Upon finding the correct facet, the next step is to determine whetherthe upwind point is still further upstream, by comparing k at theintercept with the time-step Δt. If further upwinding is required,subroutine contigf is called to find the element neighboring the newintercept facet and the logical opface is set to true, ready to callsubroutine upwelm2.F. If, however, the upwind point is within thecurrent element, the local coordinates can be determined by using Δtrather than k in equation (4), in which case opface is left false.

The upwelm2.F subroutine may be described as follows, with reference toFIG. 9. First, the local velocity, v, is computed at the currentintercept x₁, by interpolation on the facet L. Because some time hasalready been used in reaching L from the original node A, there is onlysome time Δt_(rem) remaining of the original time-step Δt. There arethree candidate facets which are checked by testing the intercept of thetrajectory with each of their planes. This can be done by settingL_(i)=0 for each node of the current facet, nodes a, b, and c in FIG. 9.

A position x is on the trajectory if:x=x ₁ −kv  (21)The local coordinates of points on the trajectory are therefore:L=J ⁻¹·(x ₁ −kv−x ₄)  (22)

The equations below are used for each node a, b, c in turn. Forsimplicity, assume the following cases are for testing the facet solelyopposite node a. In a first case, suppose node a is the fourth node inthe element node order of the tetrahedron shown in FIG. F. Set L₄=0 tofind the facet plane intercept. This implies equation (15), which isexpanded using equation (22) to give: $\begin{matrix}{{\sum\limits_{i = 1}^{3}{\sum\limits_{j = 1}^{3}{J_{ij}^{- 1}\left( {x_{I}^{j} - {kv}_{j} - x_{4}^{j}} \right)}}} = 1} & (23)\end{matrix}$

Solving this for k yields: $\begin{matrix}{k = \frac{{\sum\limits_{i = 1}^{3}{\sum\limits_{j = 1}^{3}{J_{ij}^{- 1}\left( {x_{I}^{j} - x_{4}^{j}} \right)}}} - 1}{\sum\limits_{i = 1}^{3}{\sum\limits_{j = 1}^{3}{J_{ij}^{- 1}v_{j}}}}} & (24)\end{matrix}$

In a second case, suppose a is not the fourth node, but some node m,where m is any of 1, 2, or 3. Then the intercept of the trajectory withthe plane of the facet opposite node a occurs where L_(m)=0.Substituting this in equation (22), yields: $\begin{matrix}{{\sum\limits_{j = 1}^{3}{J_{mj}^{- 1}\left( {x_{I}^{j} - {kv}_{j} - x_{4}^{j}} \right)}} = 0} & (25)\end{matrix}$

Solving this for k yields: $\begin{matrix}{k = \frac{\sum\limits_{j = 1}^{3}{J_{mj}^{- 1}\left( {x_{I}^{j} - x_{4}^{j}} \right)}}{\sum\limits_{j = 1}^{3}{J_{mj}^{- 1}v_{j}}}} & (26)\end{matrix}$

To determine if this facet is valid, the same checks as described forupwelm1.F above can be performed. Upon finding the correct facet, thenext step is to determine whether the upwind point is still furtherupstream by comparing k at the intercept with Δt_(rem). If furtherupwinding is required, subroutine contig.f can be called to find theelement neighboring the new intercept facet and the above procedure canbe repeated. In this case, Δt_(rem) is decremented appropriately. If,however, the upwind point is within the current element, the localcoordinates can be determined by using Δt_(rem), rather than k inequation (12), and upwelm2.F can then be exited.

Once the upstream point which will flow into a target node in a giventime increment has been identified in accordance with the methodsdescribed above, the temperature at that upstream point needs to bedetermined. One accurate method of determining the temperature of theupstream point employs a one dimensional analytic temperature function,as depicted in FIGS. 10A and 10B. As before, triangular elements areshown for ease of visualization, but the concept applies to threedimensional elements, without loss of generality.

In general, for each node, the respective upstream point that will“flow” or be convected to that node in the required time increment isdetermined. This may be a point in the element of which the node forms apart, or it may be a point within another element some distance away.The upstream point may or may not coincide with an upstream node. Ifnot, interpolation is required to estimate the temperature at theupstream point. According to one method, two steps may be used.

In the first step, the distance to the cavity wall of the upstream pointis estimated. While in this embodiment the wall is used as the referencesurface, in other embodiments other datums could be used, such as thesurface of the frozen layer, locally. In determining the distance to thewall, a variety of methods could be used. Since, in general, theupstream point will not be coincident with a node, its distance to thewall can be determined by linear interpolation of the distances to thewall of all the nodes in the element that enclose the upstream point.This interpolation makes use of the element shape function.

In the second step, the one dimensional analytic function at each of thenodes in the element is used to provide estimates of the temperature atthe upstream point by substituting the distance to the wall of theupstream point. The final temperature at the upstream point is obtainedby interpolation of the four node estimates according to the elementshape function. This then, is the value of temperature which, in a giventime increment, is convected to the node under consideration.

In another embodiment, the distance to the wall and the one dimensionalanalytic function describing the temperature distribution of theupstream point is taken as the nearest node in the upstream element.

While there are specific cases to be dealt with, for example, if theupstream point is the injection location, these cases do not detractfrom the fundamental idea or use of the one dimensional analyticfunction to estimate temperature.

Still another innovative methodology is related to anisotropic meshrefinement, including a distance to the mold cavity wall computationalgorithm. Anisotropic mesh refinement is a powerful means of increasingresolution in a mold cavity without an unnecessary increase incomputational burden through over refinement. The embodied refinementprocess breaks down existing tetrahedral elements into smallertetrahedra. However, the process is only applied where an insufficientnumber of layers is present. If sufficient layers are present to capturewith accuracy the solution of variables, then no refinement takes place.The algorithm runs as follows, according to one embodiment of theinvention.

An existing tetrahedral mesh is labeled with layer numbers. First, labelall nodes on the cavity wall as layer zero. Then, label all nodesadjacent to layer zero as layer one. Next, label all nodes adjacent tolayer one as layer two, and repeat with successively increasing layernumbers until all nodes are assigned a layer number. The purpose is toestablish how many layers exist in each region of the mold cavity and torefine to a required number of layers in only those regions that requireit.

In each region for which the number of layers is insufficient, thetetrahedral elements are split. The splitting process results inelements of greater aspect ratio than the original tetrahedron.Consequently, the refined mesh becomes anisotropic in that the aspectratio departs significantly from unity. Since the splitting processoccurs between successive layers, the refinement is biased towardsincreasing the density of elements and nodes in the thickness direction.This is precisely what is required for increasing accuracy insimulations since, as described earlier, the field variables vary mostsignificantly through the part thickness. FIGS. 11 through 17 showpossible ways in which a tetrahedron can be split.

More specifically, FIG. 11 is a schematic representation of tetrahedralelement face splitting combinations. All of the element splittingtemplates in the following drawings produce one of the three patterns onthe element faces depicted in FIG. 11. Either one, two, or three edgescan be refined. If this were not true, then mesh face continuity duringrefinement would not be possible. FIG. 12 is a schematic representationof a template for splitting a tetrahedral element on one edge. FIG. 13Ais a schematic representation of a template for splitting a tetrahedralelement on two adjacent edges, which entails a first split into twoelements, one of which is a four-sided base pyramid, then a subsequentsplit of the pyramid into to tetrahedral elements in accordance with thetemplate depicted in FIG. 13B. As a general rule when splitting afour-sided base pyramid, choose the diagonal of the base to split onwhich contains the node with the lowest number of the four base nodes.FIG. 13C is a template for splitting a tetrahedral element on twoopposite edges. The first split yields two tetrahedral elements whicheach require a split on one edge, using the appropriate template.

FIG. 14A is a schematic representation of a template for splitting atetrahedral element on three edges with a common shared face and FIG.14B is a schematic representation for splitting a tetrahedral element onthree edges with a common node. The latter results in a tetrahedralelement and a triangular prism element, which can be split in accordancewith the template depicted in FIG. 14C. Note that node A will be thenode with the lowest node number of the six nodes present. This ensuresthat the two quadrilaterals split by the first step are split with thediagonal including the lowest node. Also, the first step creates atetrahedral element and a quadrilateral-based pyramid, which can besubsequently split. FIG. 14D is a template for splitting a tetrahedralelement on three edges in series, of which two possible combinationsresult, depending on which nodes have the lowest node number.

FIG. 15A is a template for splitting a tetrahedral element on fouropposite edges. A first symmetrical split results in two triangularprisms, which in turn are split as depicted in FIG. 14C. FIG. 15B is atemplate for splitting a tetrahedral element on four adjacent edges. Thedivision produces two tetrahedral elements and two quadrilateral-basepyramids, which can be split using the appropriate template and rulesdiscussed hereinabove.

FIG. 16 is a template for splitting a tetrahedral element on five edges,which yields two tetrahedral elements, and a triangular prism and aquadrilateral-base pyramid, both of which can be split using theappropriate templates discussed hereinabove. Lastly, FIG. 17 is atemplate for splitting a tetrahedral element on six edges. Cutting atetrahedral element off of each corner yields a central body with eighttriangular faces. Splitting the central body in half on a quadrilateralplane yields two quadrilateral pyramids, which can be further splitaccording to the template discussed hereinabove.

Referring once again to the layering, it has been determined that anoptimal number of layers is six, and an optimal aspect ratio is fifteen.Values outside this range can be used, but when the number of layersdecreases, and the maximum aspect ratio increases, accuracy is generallyreduced. Conversely, when the number of layers increases and the aspectratio decreases, the number of elements grows quickly, and thecomputational cost becomes tends to become prohibitive.

In order to compute the distance to the wall for each internal node,nodes are numbered according to their layer and a tree constructed as anaid to finding the nearest wall location. The algorithm proceeds asfollows in accordance with one embodiment of the invention.

First, label all nodes on the cavity wall as layer zero. Then label allnodes adjacent to layer zero as layer one. Next, label all nodesadjacent to layer one as layer two, and repeat with successivelyincreasing layer numbers until all nodes are assigned a layer number.Thereafter, for each node, construct a tree of node associations for allconnected nodes that have a lower layer number down to layer numberzero, that is wall nodes. Then, project the target node on all facesformed by the layer zero nodes in the tree. The distance to the wall isrecorded as the shortest distance in the set of possible distances. Thismethod has the advantage of narrowing the possible search space forshortest distance. In the event that the node cannot be projected onto aface formed by nodes with layer number zero in its tree, then thedistance is taken as the distance to the nearest layer zero node.

This algorithm is presented diagramatically in FIGS. 18A and 18B. InFIG. 18A, the node layer numbering and the connectivity tree limits thesearch space for projection of the target node to faces “A” and “B,”which is a much smaller set than all possible faces in a cavity and istherefore highly efficient. In the case depicted in FIG. 18B, where thetarget node cannot be projected onto either face “A” or “B” in itsconnectivity tree, the distance, d, is taken to the nearest wall node.As discussed hereinabove, significant improvements in simulationpredictive accuracy have been achieved not only by consideration andapplication of conservation of energy principles, but also by addressingthe forming frozen or solidifying layer proximate the mold cavity wall.Limitations of conventional three dimensional injection molding modelingtechniques often relate to the large temperature gradients, movement ofthe solid/liquid interface boundary, large material property variationsduring solidification, and solidification of melt material duringmolding.

Loosely grouped, about five innovative methodologies related toadvantageous treatment of the solidified material have been identifiedwhich offer substantial improvements in predictive accuracy and/orcomputation time over conventional simulation methods. The first relatesto the determination of the position of the solid/liquid interface. Thenext is the formulation of elements with a linear variation of materialproperties throughout. The third relates to determination of aneffective viscosity function in elements containing the solid/liquidinterface. The next entails elimination of frozen nodes and elementsfrom the solution domain. And the fifth addresses computation of theeffective pressure in regions that have solidified.

Looking first to the determination of the position of the solid/liquidinterface, in an element through which a solid/liquid interface exists,that is to say in an element which contains both a solid region and aliquid region, material properties typically vary by several orders ofmagnitude across the interface boundary. In particular, the viscosity inthe solid region (if it is appropriate to use the term viscosity at all)is practically infinite, whilst in the liquid region, viscosity isfinite and relatively small. Accordingly, simple conventional elementassembly procedures used in the finite element method will result in aneffective overestimate of the elemental resistance to flow, because theelement stiffness will be biased heavily towards the high viscosityvalue. This is true regardless of how much of the element is actuallyfrozen. The net result is overprediction of pressure drop during flow.

In accordance with one embodiment of the invention, this problem can beresolved by making use of the one dimensional analytic function fortemperature to more accurately estimate the position of the solid/liquidinterface, and then defining an effective nodal viscosity related to thequantity of frozen material in an element. The position of the interfaceis chosen by examining the intersection of the temperature profile andthe solidification temperature of the polymer material. Thesolidification temperature, T_(s), defines the transition between solidand liquid states.

The process is represented schematically in FIG. 19. Again, for ease ofdepiction, two dimensional elements are used, but the innovative conceptcan be applied, without loss of generality, to three dimensionalelements. The drawing shows the determination of the solid/liquidinterface within an element using both the solidification temperature,T_(s), and the one dimensional analytic functions for temperaturedistribution at each of the liquid nodes.

The next innovative methodology concerns the formulation of elementswith a linear variation of material properties throughout each element.Much of common finite element practice results in a discretization ofthe domain, such that the field variables are computed at the nodepoints, and the material properties required for analysis are constantwithin an element. While this technique may provide acceptable resultsin certain less rigorous applications, according to the invention, inapplying the finite element method to fluid flow in injection molding,the equations are formulated such that both field variables (such aspressure, temperature, flow conductance and velocity) and materialproperties (such as viscosity) are defined at the nodes and are, unlessotherwise stated, interpolated within an element according to theelement shape function. Discretization in this way provides additionalaccuracy of computation, and has been found to be particularlyadvantageous for higher order, advanced interpolations.

For example, the equation for flow conductance given by Nakano, whichmakes use of the material property viscosity, can be re-written as:$\begin{matrix}{{\frac{\partial}{\partial x_{k}}\left( {\eta\frac{\partial\kappa}{\partial x_{k}}} \right)} = {- 1}} & (27)\end{matrix}$where the subscript k is an index 1, 2, or 3 denoting each coordinatedirection x, y, and z respectively.

Applying Galerkin Weighting integration over the element gives:$\begin{matrix}{{\int_{\Omega_{e}}{\frac{\partial}{\partial x_{k}}\left( {\eta^{e}\frac{\partial\kappa^{e}}{\partial x_{k}}} \right)N_{i}{\mathbb{d}\Omega_{e}}}} = {- {\int_{\Omega_{e}}{N_{i}{\mathbb{d}\Omega_{e}}}}}} & (28)\end{matrix}$where the elemental flow conductance and viscosity respectively aregiven by: $\begin{matrix}{\kappa^{e} = {\sum\limits_{j = 1}^{4}{\kappa_{j}N_{j}}}} & (29) \\{\eta^{e} = {\sum\limits_{j = 1}^{4}{\eta_{j}N_{j}}}} & (30)\end{matrix}$and where Ω_(e) refers to the volume of an element, N_(i) is the elementshape function at each node in the element i, and the superscript e isan element index.

Applying the divergence theorem gives: $\begin{matrix}{{{\int_{\Gamma_{e}}{\eta^{e}\frac{\partial\kappa^{e}}{\partial x_{k}}n_{k}N_{i}{\mathbb{d}\Gamma_{e}}}} - {\int_{\Omega_{e}}{\eta^{e}\frac{\partial\kappa^{e}}{\partial x_{k}}\frac{\partial N_{i}}{\partial x_{k}}{\mathbb{d}\Omega_{e}}}}} = {- {\int_{\Omega_{e}}{N_{i}{\mathbb{d}\Omega_{e}}}}}} & (31)\end{matrix}$where n_(k) is the normal vector at the element boundary Γ.

The first term is assumed to vanish during element assembly at internalelement faces due to the equal and opposite contribution from adjacentelements. This is true, since element viscosity, which is based on aninterpolation of nodal values, will be matched along the boundarybetween two adjacent elements. Further, the gradients of flowconductance are assumed to be equal on each side of the elementboundary. Only the normal vector has an opposite sign in adjacentelements, which results in the contributions from adjacent elementsnegating to zero. At the cavity boundary in the filled region, theboundary integral is not needed when a known value boundary condition isapplied (κ=0). However in the case where wall slip is required, so toois the boundary integral term. In either case the boundary integral termneeds to be evaluated at the free surface boundary and at the injectionboundary.

Equations (29) and (30), which demonstrate the linear variation ofproperties within an element, can be substituted into equation (31)giving: $\begin{matrix}{{{- {\int_{\Gamma_{e}}{\eta_{l}N_{l}\frac{\partial N_{j}}{\partial x_{k}}n_{k}N_{i}\kappa_{j}{\mathbb{d}\Gamma_{e}}}}} + {\int_{\Omega_{e}}{\eta_{l}N_{l}\frac{\partial N_{j}}{\partial x_{k}}\frac{\partial N_{i}}{\partial x_{k}}\kappa_{j}{\mathbb{d}\Omega_{e}}}}} = {\int_{\Omega_{e}}{N_{i}{\mathbb{d}\Omega_{e}}}}} & (32)\end{matrix}$It is known that the exact integral for a tetrahedron with naturalco-ordinate shape functions is: $\begin{matrix}{{\int_{\Omega_{e}}{N_{i}^{\alpha}N_{j}^{\beta}N_{k}^{\gamma}N_{l}^{\delta}{\mathbb{d}\Omega_{e}}}} = {\frac{{\alpha!}{\beta!}{\gamma!}{\delta!}}{\left( {\alpha + \beta + \gamma + \delta + 3} \right)!}6V}} & (33)\end{matrix}$where:

-   N_(i) is the shape function,-   V is the volume of the element, and-   α, β, and γ, are integer exponents

Using this integration relationship (33) with equation (32) gives:$\begin{matrix}{\begin{matrix}\left( {{\sum\limits_{f = 1}^{4}{\left( {\delta_{fi} - 1} \right)\left\lbrack {\frac{A{\sum\limits_{l = 1}^{3}{\left( {1 + \delta_{il}} \right)\eta_{l}}}}{12}\frac{\partial N_{j}}{\partial x_{k}}n_{k}} \right\rbrack}} +} \right. \\{\left. \left\lbrack {\frac{\sum{\eta_{l}V}}{4}\frac{\partial N_{j}}{\partial x_{k}}\frac{\partial N_{i}}{\partial x_{k}}} \right\rbrack \right)\kappa_{j}}\end{matrix} = \frac{V}{4}} & (34)\end{matrix}$where the summation of viscosities in the boundary integral is only overthe three nodes which form each face. For completeness, the integral isshown here on every face, f. The face numbering matches the number ofthe opposite vertex.

So it is shown how equation for flow conductance, which requires thematerial property viscosity, can be formulated as a finite elementsystem with viscosity varying throughout each element according to theelement shape function.

Still another area of innovation concerns the determination of aneffective viscosity function in elements containing the solid/liquidinterface. The most common partially-frozen elements occur adjacent tothe mold wall, where they contain both “cold” wall nodes and “hot”internal cavity nodes. A linear variation of viscosity, which relies onthe assumption of a linear variation of viscosity through elements, andarises from the use of elements with linear shape functions, usuallyover-estimates the effective viscosity of such elements, creatingsignificant error in many cases, the most common and most serious beingpartially-frozen elements. See, for example, FIG. 20, which compares alinear viscosity profile with a more realistic, highly non-linearprofile. This linear assumption leads to artificially high flowresistance in the simulation, reflected in injection pressures which canbe more than 100% greater than reality. The approximation of linearviscosity variation is thus very poor in partially-frozen elements.While the error asymptotes to zero for a sufficiently fine mesh, thissolution is not generally computationally practical or efficient.Accordingly, absent implementation of an effective viscosity functiontechnique, simulation accuracy can suffer greatly.

According to one embodiment of the invention, greatly improvedsimulation accuracy of viscosity in partially-frozen elements entailscalculating an effective viscosity. More specifically, the effectiveviscosity at each frozen node in an element that contains thesolid/liquid interface can be computed and used locally during theelement assembly procedure. As such, these nodal viscosities aretransient, and do not replace the global values computed through aknowledge of the temperature and shear rate state at each node and theviscosity material function.

By labeling nodes above and below a predetermined solidificationtemperature, T_(s), as “hot” and “cold” respectively, partially-frozenelements are defined as those which are not completely “hot” (i.e. allnodes designated “hot”) or completely “cold” (i.e. all nodes designated“cold”). Cold nodes may be assigned an arbitrary large viscosity, thatis, a viscosity value sufficiently large to ensure that no flow occursin frozen portions. The solid state viscosity may be typically on theorder of 10⁶ MPas. This technique gives a more realistic effectiveviscosity in partially-frozen elements. The viscosities of completelyfrozen or completely hot elements need not be modified.

Before addressing the specifics of one embodiment of the effectiveviscosity method, it may be helpful to provide an overview of some ofthe main features thereof. The modifications to viscosities can be madeto any method for solving non-isothermal flow, including the method byNakano, and more conventional Stokes and Navier-Stokes solutions. Whenshear heating is enabled, the appropriate viscosity modifications can becarried through related routines.

Additionally, the modifications are made during element assembly. Foreach element, the four local node viscosities are modified. Global nodeviscosities are not altered. Viscosities in elements whose nodes are allcold or all hot are not modified. The standard approximation of linearviscosity variation is still used for these elements.

The volume fraction of freeze in the element is used to calculate avolume fraction of freeze within each nodal volume. In this particularcontext, “freeze” is the portion of the volume below T_(s) and nodalvolume refers only to that part within the current element. The meanviscosities, vis_(h) and vis_(c) of the hot and cold nodes in theelement, respectively, are calculated and each nodal volume may containa frozen and a molten portion, regardless of whether the node itself ishot or cold. The viscosity within the frozen portion is the nodalviscosity if the node is cold, and vis_(c) if the node is hot.Similarly, the viscosity within the molten portion is the nodalviscosity if the node is hot, and vis_(h) if the node is cold. Lastly,the effective viscosity of the entire nodal volume is then determinedand applied locally to the node. The expression for effective viscositywas derived by determining the parallel flow resistance of the frozenand molten portions of the nodal volume, based on the assumption of flowparallel to the local freeze surface.

Table 1 lists the terms used in this methodology, for ease of reference.TABLE 1 Variable Definition T_(s) Solidification temperature nc Numberof nodes on element below T_(s) nh Number of nodes on element aboveT_(s) tfroz Volume fraction of freeze in element (freeze meaning belowT_(s)) frzc Volume fraction of freeze in nodal volume of each cold nodefrzh Volume fraction of freeze in nodal volume of each hot node vis_(i)Viscosity of element nodes (i = 1 to 4) — global values etet_(i)Viscosity of element nodes (I = 1 to 4) — local values used for elementassembly vis_(c) Mean of global viscosities of nc cold nodes vis_(h)Mean of global viscosities of nh hot nodes

For each element, the following steps are followed to determine thelocal node viscosities, etet_(i). Derivation of the equations isdescribed hereinbelow. According to this embodiment of the method, firstdetermine nc and nh. If nc=0 or nh=0, this is not a partially-frozenelement. Set etet_(i)=vis_(i) and proceed to the next element. Next,calculate vis_(c) and vis_(h). Then, calculate frzc and frzh inaccordance with the following equations: $\begin{matrix}{{frzc} = {{Min}\left( {\frac{4 \cdot {tfroz}}{nc},1} \right)}} & (35) \\{{frzh} = \frac{{4 \cdot {tfroz}} - {{nc} \cdot {frzc}}}{nh}} & (36)\end{matrix}$

Finally, loop over the four local nodes i and calculate the viscosityetet_(i) of each. Use equation (37) if node i is cold, and equation (38)if node i is hot. $\begin{matrix}{{{etet}_{i} = \frac{{vis}_{i}{vis}_{h}}{{\left( {1 - w} \right){vis}_{i}} + {w \cdot {vis}_{h}}}},{w = \frac{{frzc}^{1/3}}{{frzc}^{1/3} + \left( {1 - {frzc}} \right)^{1/3}}}} & (37) \\{{{etet}_{i} = \frac{{vis}_{i}{vis}_{c}}{{\left( {1 - w} \right){vis}_{c}} + {w \cdot {vis}_{i}}}},{w = \frac{{frzh}^{1/3}}{{frzh}^{1/3} + \left( {1 - {frzh}} \right)^{1/3}}}} & (38)\end{matrix}$

The quantities frzc and frzh are derived as follows. First, let thecurrent element volume be V, so that the four nodal volumes are eachV/4. The volume of freeze in the element is V*tfroz. By approximatingthat the freeze volume is distributed evenly between the cold nodes, nc,there is a volume V*tfroz/nc of freeze in each nodal volume. Since thisfreeze volume cannot be greater than the nodal volume, the actual freezevolume is the minimum of V*tfroz/nc and V/4, ie. Min(V*tfroz/nc, V/4).Converting this to a volume fraction yields equation (35), namely:$\begin{matrix}{{frzc} = {{Min}\left( {\frac{4 \cdot {tfroz}}{nc},1} \right)}} & (35)\end{matrix}$

Any freeze “left over” or remaining after distribution to the cold nodesis distributed evenly to the hot nodes, nh. The freeze volume assignedto cold nodes is V*nc*frzc/4, so the remaining freeze volume isV*tfroz−V*nc*frzc/4. Dividing this between the hot nodes, nh, andconverting to a volume fraction yields equation (36), namely:$\begin{matrix}{{frzh} = \frac{{4 \cdot {trfoz}} - {{nc} \cdot {frzc}}}{nh}} & (36)\end{matrix}$

The local node viscosities are derived as follows. Let the nodal volumebe Vt. Assuming the volume is partially frozen, the total volume isdivided into a cold volume, V_(c), and a hot volume, V_(h), withrespective viscosities η_(c) and η_(h). Obviously, V_(t)=V_(c)+V_(h).Also, the local flow should be parallel to the local freeze surface inthe nodal volume. FIG. 21 depicts this situation schematically. Whilethe figure is in two dimensions, the drawing is representative of threedimensions.

The next step is to determine the effective viscosity of the entirenodal volume V_(t) in the flow direction. The effective viscosity normalto the flow does not affect the flow solution. Call this effectiveviscosity η_(t). Exact calculation of η_(t) tends to be impractical forodd-shaped volumes in the model, requiring complex volume integrations.Instead, a simplified approach can be used which gives the correctbehavior as each nodal volume freezes.

Now, the flow resistance R of some volume with uniform viscosity η isproportional to ηL/A, where L is the path length through the volume inthe flow direction and A is the cross-sectional area normal to the flowdirection. The flow resistance of the nodal volume is given by1/R_(t)=1/R_(c)+1/R_(h). Calculating L/A for each volume is generallynot feasible. Instead, note that L/A has the dimensions of 1/length, soapproximate that 1/A varies as volume^(−1/3). Therefore, the followingrelationship applies: $\begin{matrix}{\frac{V_{t}^{1/3}}{\eta_{t}} = {\frac{V_{c}^{1/3}}{\eta_{c}} + \frac{V_{h}^{1/3}}{\eta_{h}}}} & (39)\end{matrix}$from which can be obtained: $\begin{matrix}{\eta_{t} = \frac{\eta_{c}\eta_{h}V_{t}^{1/3}}{{\eta_{h}V_{c}^{1/3}} + {\eta_{c}V_{h}^{1/3}}}} & (40)\end{matrix}$

Let the frozen fraction in the nodal volume be F, where F is frzc orfrzh, depending on whether the node is cold or hot. Then V_(c)=FV_(t)and V_(h)=(1−F)V_(t). Substituting into equation (40) yields:$\begin{matrix}{\eta_{t} = \frac{\eta_{c}\eta_{h}}{{\eta_{h}F^{1/3}} + {\eta_{c}\left( {1 - F} \right)}^{1/3}}} & (41)\end{matrix}$

Equation (41) may be normalized. For example, if η_(c)=η_(h), equation(41) should reduce to η_(t)=η_(c). Accordingly modify equation (41) to:$\begin{matrix}{\eta_{t} = \frac{\eta_{c}{\eta_{h}\left\lbrack {F^{1/3} + \left( {1 - F} \right)^{1/3}} \right\rbrack}}{{\eta_{h}F^{1/3}} + {\eta_{c}\left( {1 - F} \right)}^{1/3}}} & (42)\end{matrix}$

The final question is what viscosity to use in the two portions of thenodal volume. If the node i associated with the nodal volume is cold,use η_(c)=vis_(i) and approximate that η_(h)=vis_(h). Similarly, if nodei is hot, use η_(h)=vis_(i) and approximate that η_(c)=vis_(c). Thefinal equation for etet_(i) is then equation (37) or equation (38)according to whether the node i is cold or hot.

The behavior of equation (42) is correct in limiting cases. For example,if F=0, η_(t)=η_(h) as required. Similarly, if F=1, η_(t)=η_(c) asrequired. Also, the variation of η_(t) with F is smooth and monotonic.If η_(c)=η_(h), η_(t)=η_(c) as required. Lastly, when η_(c)>>η_(h),η_(t)≈η_(h)[1+(F/1−F^(1/3)]. This has the correct form, as lit increaseswith F from a value of η_(h) at F=0.

Two simpler expressions in place of equation (42) can be used, ifdesired; however, they may result in less accurate predictive behavior.One expression is as follows:η_(t)=η_(c) F+η _(h)(1−F)  (43)

One issue with equation (43) is that in the typical case thatη_(c)>>η_(h), the equation yields η_(t)≈η_(c)F. Effectively, since η_(c)is very large, this equation predicts there is no flow through the node,even when most of the node volume is molten. Correspondingly, thepredicted injection pressures tend to be larger than expected.

The second expression is as follows: $\begin{matrix}{\eta_{t} = \frac{\eta_{c}\eta_{h}}{{\eta_{h}F} + {\eta_{c}\left( {1 - F} \right)}}} & (44)\end{matrix}$

This equation does give the correct behavior in limiting cases, likeequation (42). Namely if F=0, η_(t)=η_(h) as required. If F=1,η_(t)=η_(c) as required. The variation of η_(t) with F is smooth andmonotonic. If η_(c)=η_(h), η_(t)=η_(c) as required. Lastly, whenη_(c)>>η_(h), it η_(h)[1+(F/1−F)]. This has the correct form, as itincreases with F from a value of η_(h) at F=0.

There are two issues, however, with equation (44). Firstly, equation(44) is derived by assuming flow resistance proportional to η/V, ratherthan η/V^(1/3) as for equation (42). The underlying assumption istherefore dimensionally incorrect. Secondly, in practice, equation (42)has been demonstrated to predict injection pressures closer to thecorrect values. For a tube model test case, with an expected injectionpressure of 0.62 MPa, 6-layer Darcy models gave pressures of 0.47 MPausing equation (44), compared to 0.51 MPa using equation (42).Sixteen-layer models gave pressures of 0.56 MPa and 0.55 MParespectively. A linear viscosity formulation gave pressures of 1.59 MPaand 0.63 MPa with 6 and 16 layer refinement.

Testing of the effective viscosity methodology has shown that the valueof tfroz is relatively accurate. Further, while the present methodemploys equal distribution of tfroz to the “cold” nodes, an alternativeapproximation would be to weight the distribution by the element nodaltemperatures. Yet another would be to directly calculate frzc byinterpolating the T_(s) surface through each nodal volume andcalculating volumes. A similar technique could be employed for the equaldistribution of remaining freeze to the “hot” nodes, although thesealternatives could consume additional computing resources.

While the use of F^(1/3) and (1−F)^(1/3) weightings for cold and hotregion viscosities are based upon empirical weightings, with additionalqualitative theoretical support, it is contemplated that the weightingscould be based on integration over partial nodal volumes. Additionally,the use of mean hot node viscosity for the hot region of a cold node maytend to underestimate the hot region viscosity, since viscosityincreases as the freeze surface is approached. This, in turn, may resultin a slight underestimation of injection pressure, as noted hereinabove.The use of mean cold node viscosity for the cold region of a hot node,however, has been determined to produce negligible error, since vis_(c)simply needs to be sufficiently large to eliminate flow through thefrozen portion.

Now, with respect to the technique of elimination of frozen nodes andelements from the solution domain, at some stage during filling and/orpacking, whole elements are rendered substantially immobile by virtue ofthe fact that the temperature of all nodes in the element has fallenbelow the solidification temperature, T_(s). The question is how totreat these elements within the context of fluid flow analysis. Oneapproach is to continue to treat these elements as a fluid byconsidering the solid region as a zone of very high viscosity. The“solid” viscosity can be obtained by extrapolating the fluid viscosityfunction, or by assigning an arbitrarily large value. This, however,this can lead to ill-conditioning of the system of equations.Convergence can be retarded, or stopped completely, and thecomputational cost increased thereby.

In accordance with one embodiment of the invention, the solution domainfor the flow conductance and pressure is reduced to include only thoseelements which have some nodes above the solidification temperature. Byremoving completely frozen elements from the solution domain, ease ofconvergence for an essentially fluid problem is achieved, while a savingin computation time is also gained. The useful boundary conditions forno-slip flow are applied at the boundary with the frozen region, ratherthan at the mold wall. In the early phases of the molding process, thesetwo boundaries are coincident. These boundary conditions are applied bysetting the flow conductance function to be zero at the boundary. Notethat the full cavity domain is typically used for the temperature fieldcalculations.

Still another innovative methodology entails computation of theeffective pressure in regions that have solidified. Having chosen tocharacterize the material in the mold cavity as either solid or liquid,it is still generally desirable to compute the pressure in the solid. Inreality, the frozen polymer still experiences a pressure influence,since compressive stress and small deflections can be transmitted by thesolid polymer from the neighboring molten polymer over small distances.In order to assign a pressure to the frozen nodes, it is generally notsufficient to propagate pressure from molten nodes, because this methoddoes not allow for the process of freeze-off, where some sections of thefilled cavity can still be molten and subjected to injection pressure,while at the same time other regions are frozen and have no appliedpressure because they are at such a distance from the molten regionsthat no meaningful pressure transmission is effected.

Also, it is desirable to know the pressure in order to establish thestate of the material. The pressure is used along with a knowledge oftemperature to establish material density and, consequently, mass ateach node, solid or liquid, within the cavity.

Further, it is reasonably common for manufacturers to flush mountpressure sensors within actual mold cavities in order to provide someinformation by which they can control molding machines and processes,and attempt to establish parameters related to molded part quality.Consideration of the actual physics of the process leads one to concludethat a precise simulation likely would require a very complex andcomputationally expensive multi-physics simulation, combining both fluiddynamics and solid mechanics. At this time, however, the techniques arenot well established.

One solution, in accordance with the invention, has been to develop atechnique for estimating solid pressure by projecting core pressure ontothe outer cavity frozen layer. More specifically, in both frozen andmolten regions, the core pressure is projected onto any outer nodeswhich are frozen. By way of example, the core pressure may be determinedfrom the normal fluid finite element method calculation, or it may befrom the pressure decay of the frozen solid material, which progressesas the temperature of the polymer decreases due to cooling. By a processof geometry decomposition, all nodes which are not in the core of a partare assigned a core node, whose pressure they will mimic once frozen.The only nodes for which pressure decay during cooling is actuallycalculated then, are those which are given the core attribute. Eachouter node may be assigned the core node which is closest to it to beits core node. However, the algorithm for determining which nodes arecore nodes actually relies on the outer nodes. That is, every core nodehas at least one outer node which depends on it; otherwise, the corenode would not have been given the core node attribute.

The algorithm for generating the core node attribute is as follows.First, label all nodes on the cavity wall as layer zero. Then, label allnodes adjacent to layer zero as layer one. Next, label all nodesadjacent to layer one as layer two, and repeat with successivelyincreasing layer numbers until all nodes are assigned a layer number.Thereafter, calculate the distance to the nearest cavity wall for everynode. Also, record which is the closest cavity wall node for eachinternal node. All nodes whose layer number is greater than or equal tohalf the guaranteed number of layers through the thickness are marked aspotential core nodes.

The guaranteed number of layers through the thickness is a result of themesh refinement process described earlier. Now, continue by startingwith each node on the cavity wall, that is, with layer number zero.Traverse the nodes connected to each in turn, moving to successivelyhigher layer numbers, until such time as a potential core node isreached. In general, there will be several potential core nodes for eachcavity wall node. The node that is selected as the core node associatedwith the wall node is that which maximizes the ratio d_(cw)/d_(cn),where:

-   -   d_(cw) is distance of core node to cavity wall, and    -   d_(cn) is distance of core node to wall node.        Internal nodes that are not core nodes are assigned a core node        that is the same as the closest cavity wall node, determined        during the previous layer numbering steps.

This algorithm has the feature that when a thin and thick region areadjacent to each other, the frozen nodes in the thin region will have acore pressure from the thick region assigned to them. This has beenfound to mimic the real physical process well. For example, where anarrow gate is attached to a thicker runner feed system, after the gatehas frozen, the pressure reported on it mimics that of the thicker feedsystem. This represents the transmission of compressive stresses throughthe solid polymer.

FIG. 22 shows layer numbering for core, wall, and internal nodes and thedefinition of the distances used in the ratio to establish which corenode is assigned to each wall node.

FIGS. 23A through 23C shows the relationship between core nodes, wallnodes and internal nodes, and their change in pressure with increasingtime and freezing.

Having described the various innovative methodologies and algorithms indetail, attention can now be turned to FIGS. 24 and 25 for anunderstanding of how the improvements are implemented in both thefilling phase and packing phase, in accordance with one embodiment ofthe invention. FIG. 24 is a schematic representation of the fillingphase analysis summarizing certain process substeps of step 50 from FIG.2. Recall, that according to the top level flowchart of FIG. 2, once thethree dimensional solid model has been discretized to provide a modelsolution domain and the boundary conditions have been set, thesimulation first solves for filling phase process variables in step 50.

The filling phase analysis of step 50 begins with the assumption thatthe cavity is empty in step 110. All field variables, such as pressure,temperature, and velocity, are initialized in step 120. Thereafter,conservation of mass and conservation of momentum equations are appliedto solve for fluidity, κ, in step 130 and pressure in step 140 andcalculate velocity in step 150, for at least a portion of the solutiondomain. In the case where a finite element analysis is being used, thefluidity, pressure, and velocity values are determined at each of thenodes in the mesh.

Fluidity can be solved by application of the Darcy flow equation andpressure then solved using a Laplace equation, because pressure is afunction of the fluidity. Alternatively, a Navier-Stokes technique canbe used to solve for pressure, based on conservation of momentum. Oncefluidity and pressure have been solved, velocity can be calculateddirectly. Next, based on conservation of energy principles, temperatureis used to calculate viscosity in step 160. Then, in step 170, the modelchecks to see if pressure has converged, which is likely not the caseinitially. If not, steps 130 through 160 are repeated iteratively, untilpressure converges.

Once pressure has converged, the simulation incrementally advances thefluid flow front in free surface evolution step 180. The contribution ofconvective heat transfer is next considered in step 190, in which anexplicit Lagrangian temperature analysis is completed to account for thethermal energy effects of the continuing inflow of molten polymer in thecavity. Then, the remaining terms in the temperature equation are solvedin step 200, to account for conduction, viscous dissipation, and anyother thermal effects desired, such as heat of solidification inthermoplastics and heat of reaction in thermosetting materials.

Once all of the temperature effects have been quantified and integratedinto the simulation, the algorithm checks to see if the mold cavity isfull in step 210. If not, steps 130 through 200 are repeated,iteratively, until such time as the cavity is full. In general, thealgorithm may loop on the order of 100 to 200 times to simulate fillingof the mold cavity. Upon completion, the model then proceeds to step 60,depicted in the flowchart in FIG. 2, wherein the packing phase analysisis performed.

Referring now to FIG. 25, depicted is a schematic representation of thepacking phase analysis summarizing certain process substeps of step 60from FIG. 2. The packing phase analysis of step 60 begins with theinitial state of all variables resulting from the filling phase, step50. Once again, conservation of mass and conservation of momentumequations are applied to solve for fluidity, κ, in step 220 and pressurein step 230 and calculate velocity in step 240, for at least a portionof the solution domain. Once fluidity and pressure have been solved,velocity can be calculated directly. Next, based on conservation ofenergy principles, temperature is used to calculate viscosity in step250. Then, in step 260, the model checks to see if pressure hasconverged, which is likely not the case initially. If not, steps 220through 150 are repeated iteratively, until pressure converges.

Once pressure has converged, the contribution of convective heattransfer is next considered in step 270, in which the explicitLagrangian temperature analysis is completed to account for the reduced,but continuing inflow of molten polymer in the cavity. Then, theremaining terms in the temperature equation are solved in step 280, toaccount for conduction, viscous dissipation, and any other thermaleffects.

Once all of the temperature effects have been quantified and integratedinto the simulation, component properties such as density, volumetricshrinkage, mass, and frozen volume are calculated and updated in step290. Thereafter the algorithm checks to see if the prescribed pressureprofile has been completed in step 300. If not, steps 220 through 290are repeated, iteratively, until such time as the pressure profile iscomplete. Upon completion, the model then proceeds to step 70, depictedin the flowchart in FIG. 2, wherein the results of the simulation areoutput for consideration by the design engineer.

While many of the steps in the filling phase and packing phaseflowcharts are similar, the respective contributions to the predictiveaccuracy of the overall simulation vary, especially with respect toapplication of the principles of conservation of energy. For example,during the filling phase where there exists relatively high flow rates,viscous dissipation thermal effects can be significant; whereas, duringthe packing phase, where there exists relatively low flow rates, thecontributions of conduction and convection dominate the heat transferanalysis.

With regard to the numerous methodologies discussed hereinabove relatedto the principles of conservation of energy for improving predictiveaccuracy of the model, according to one embodiment of the invention, thefollowing techniques may be implemented advantageously in step 190 ofthe filling phase flowchart: the use of a one dimensional analyticfunction to describe the local temperature distribution at a node; thedefinition of the variation of a one dimensional analytic function, withtime, to account for heat convection; the description of an explicittemperature convection scheme using the one dimensional analyticfunction; and finite element mesh refinement, including a distance tothe cavity wall computation algorithm. The description of the variationof a one dimensional analytic function to account for viscous heatgeneration may be implemented in step 200. These techniques could alsobe implemented in corresponding steps in the packing phase flowchart;however, their influence in the simulation results would be to a lesserdegree, due to the different conditions.

Similarly, with respect to the innovative methodologies related totreatment of the solidified material, the following techniques may beimplemented advantageously in step 220 of the filling phase flowchart insolving for fluidity: the determination of the position of thesolid/liquid interface; the formulation of elements with a linearvariation of material properties throughout; and the determination of aneffective viscosity function in elements containing the solid/liquidinterface. The elimination of frozen nodes and elements from thesolution domain and the computation of the effective pressure in regionsthat have solidified may be implemented effectively in step 200 tofacilitate solving for pressure. These techniques could also beimplemented in corresponding steps in the filling phase flowchart;however, their influence in the simulation results would be to a lesserdegree.

While there has been described herein what are considered to beexemplary and preferred embodiments of the invention, othermodifications and alternatives of the invention will be apparent tothose skilled in the art from the teachings herein. All suchmodifications and alternatives are considered to be within the scope ofthe invention. For example, while the disclosure has been generallydirected to the modeling of thermoplastic injection molding, theteachings of the invention are applicable to injection molding ofthermosetting polymers. As mentioned above, in such simulations, thecontribution of the heat of reaction may be considered in theconservation of energy analysis, as can the heat of solidification inthe thermoplastic model. Other heat transfer and thermal energyconsiderations can be integrated into the analysis, such as the thermalcapacitance of the mold and active cooling or heating of the mold.

The general principles underlying the teachings herein may also beapplied to modeling flows of different materials simultaneously orsequentially in the same mold, as well as modeling of other materialsand fluids. For example, some or all of the principles of the inventionpresented herein could be applied to other molding processes. This couldinclude polymer processes such as extrusion, blow molding, compressionmolding, thermoforming, and injection molding variants such as gasassisted injection molding, compression injection molding, andco-injection molding. Moreover, the principles could also be applied tomanufacturing processes in other industries such as metal casting andfood processing, namely, anywhere where a fluid, or a substance whichbehaves like a fluid, flows.

Accordingly, what is desired to be secured by Letters Patent is theinvention as defined and differentiated in the following claims andequivalents thereof.

1-25. (canceled)
 26. A method for modeling injection of a fluid into amold defining a three dimensional cavity, the method comprising thesteps of: (a) providing a three dimensional computer model defining thecavity; (b) discretizing a solution domain based on the model; (c)specifying boundary conditions; (d) solving for filling phase processvariables over at least part of the solution domain to providerespective filling phase solutions therefor; and (e) solving for packingphase process variables over at least part of the solution domain usingrespective states of the process variables at termination of filling; toprovide respective packing phase solutions therefor, wherein at leastone of steps (d) and (e) comprises the substeps of: using a firstdescription of a distribution of a process variable about each of aplurality of nodes or elements within the solution domain; and using asecond description of the distribution of the process variable in atleast part of the solution domain comprising the plurality of nodes orelements, the second description comprising conservation of mass,conservation of momentum, and conservation of energy equations.
 27. Themethod according to claim 26, wherein the filling phase processvariables and packing phase process variables are selected from thegroup consisting of density, fluidity, mold cavity fill time, moldcavity packing time, pressure, deformation rate, shear stress,temperature, internal energy, velocity, velocity gradient, flow rate,viscosity, and volumetric shrinkage.
 28. The method according to claim26, further comprising the steps of: (f) determining whether at leastone of the respective filling phase solutions and packing phasesolutions are acceptable; (g) modifying at least one of the discretizedsolution domain and the boundary conditions in the event at least one ofthe respective filling phase solutions and packing phase solutions isdetermined to be unacceptable; and (h) repeating steps (d) through (g)until the respective filling phase solutions or packing phase solutionsare determined to be acceptable.
 29. The method according to claim 26,further comprising the step of: displaying in graphics format a fillingphase solution selected from the group consisting of fill time,pressure, deformation rate, shear stress, temperature, velocity, andviscosity.
 30. The method according to claim 26, further comprising thestep of: displaying in graphics format a packing phase solution selectedfrom the group consisting of density, packing time, pressure,deformation rate, temperature, velocity, viscosity, and volumetricshrinkage.
 31. A method for modeling injection of a fluid into a molddefining a three dimensional cavity, the method comprising the steps of:(a) providing a three dimensional computer model defining the cavity;(b) discretizing a solution domain based on the model; (c) specifyingboundary conditions; and (d) solving for filling phase process variablesover at least part of the solution domain to provide respective fillingphase solutions therefor, wherein step (d) comprises the substeps of:using a first description of a distribution of a process variable abouteach of a plurality of nodes or elements within the solution domain; andusing a second description of the distribution of the process variablein at least part of the solution domain comprising the plurality ofnodes or elements, the second description comprising conservation ofmass, conservation of momentum, and conservation of energy equations.32. The method according to claim 31, wherein the discretizing step (b)comprises the substep of generating a finite element mesh based on themodel by subdividing the model into a plurality of connected elementsdefined by a plurality of nodes.
 33. The method according to claim 31,wherein the boundary conditions are selected from the group consistingof fluid composition, fluid injection location, fluid injectiontemperature, fluid injection pressure, fluid injection volumetric flowrate, mold temperature, cavity dimensions, cavity configuration, andmold parting plane, and variations thereof.
 34. The method according toclaim 31, wherein the solving step (d) utilizing the conservation ofmass and conservation of momentum equations comprises the substeps of:(i) solving for fluidity over at least part of the solution domain; (ii)solving for pressure over at least part of the solution domain; and(iii) calculating velocity over at least part of the solution domain.35. The method according to claim 34, wherein the solving step (d)utilizing the conservation of energy equation comprises the substep ofcalculating viscosity over at least part of the solution domain.
 36. Themethod according to claim 35, wherein the viscosity calculating substepis based on temperature.
 37. The method according to claim 36, whereinat least one of velocity and viscosity is calculated iteratively, untilpressure converges.
 38. The method according to claim 37, wherein thesolving step (d) comprises the substep of determining free surfaceevolution of the fluid in the cavity based on velocity.
 39. The methodaccording to claim 38, wherein the solving step (d) comprises thesubstep of calculating temperature based on at least one of a convectiveheat transfer contribution, a conductive heat transfer contribution, anda viscous dissipation contribution.
 40. The method according to claim39, wherein free surface evolution is determined iteratively, until thecavity is filled.
 41. The method according to claim 31 furthercomprising the step of: (e) solving for packing phase process variablesusing conservation of mass, conservation of momentum, and conservationof energy equations for at least a portion of the solution domain basedin part on respective states of the process variables at termination offilling, to provide respective packing phase solutions therefor for atleast some of the portion of the solution domain.
 42. The methodaccording to claim 41, wherein the solving step (e) utilizing theconservation of mass and conservation of momentum equations comprisesthe substeps of: (i) solving for fluidity over at least part of thesolution domain; (ii) solving for pressure over at least part of thesolution domain; and (iii) calculating velocity over at least part ofthe solution domain.
 43. The method according to claim 41, wherein thesolving step (e) utilizing the conservation of energy equation comprisesthe substep of calculating viscosity over at least part of the solutiondomain.
 44. The method according to claim 43, wherein the viscositycalculating substep is based on temperature.
 45. The method according toclaim 44, wherein at least one of velocity and viscosity is calculatediteratively, until pressure converges.
 46. The method according to claim45, wherein the solving step (e) comprises the substep of calculatingtemperature based on at least one of a convective heat transfercontribution, a conductive heat transfer contribution, and a viscousdissipation contribution.
 47. The method according to claim 46, furthercomprising the step of: (f) calculating mass properties of a component.48. The method according to claim 47, wherein the mass properties areselected from the group consisting of component density, volumetricshrinkage, component mass, and component volume.
 49. The methodaccording to claim 47, wherein at least one of velocity, viscosity, andmass properties is calculated iteratively, until a predeterminedpressure profile is completed.
 50. The method according to claim 32,wherein the mesh generating substep comprises generating an anisotropicmesh in thick and thin zones such that mesh refinement providesincreased resolution in a thickness direction without increasingsubstantially mesh refinement in a longitudinal direction.
 51. Themethod according to claim 26, wherein at a given time step, the firstdescription describes each of the plurality of nodes or elementsindependently of the others.
 52. The method according to claim 26,wherein the value of the process variable at a first point provided bythe first description of the process variable about a first node orelement is not necessarily equal to the value of the process variable atthe first point provided by the first description of the processvariable about a second node or element.
 53. The method according toclaim 26, wherein the value of the process variable at a first pointprovided by the first description of the process variable about a nodeor element is used in the second description.
 54. The method accordingto claim 26, wherein the first description is a one dimensional analyticfunction or is a discrete function.
 55. The method according to claim26, wherein the first description describes a distribution oftemperature or internal energy about a node or element.
 56. The methodaccording to claim 55, wherein the first description is or approximatesa solution for one dimensional heat conduction in a solid.
 57. Themethod according to claim 26, wherein the first description comprises oris derived from an error function.
 58. The method according to claim 26,wherein the first description is a one dimensional description oftemperature distribution about a node or element, the descriptioncomprising an error function.
 59. The method according to claim 31,further comprising the step of: (e) determining whether the respectivesolutions are acceptable for injection of the fluid during filling ofthe mold cavity.
 60. The method according to claim 31, wherein at agiven time step the first description describes each of the plurality ofnodes or elements independently of the others.
 61. The method accordingto claim 31, wherein the value of the process variable at a first pointprovided by a first description of the process variable about a firstnode or element differs from the value of the process variable at thefirst point provided by a first description of the process variableabout a second node or element.
 62. The method according to claim 31,wherein the value of the process variable at a first point provided bythe first description of the process variable about a node or element isused in the second description.
 63. The method according to claim 31,wherein the first description is a one dimensional analytic function oris a discrete function.
 64. The method according to claim 31, whereinthe first description describes a distribution of temperature orinternal energy about a node or element.
 65. The method according toclaim 64, wherein the first description is or approximates a solutionfor one dimensional heat conduction in a solid.
 66. The method accordingto claim 31, wherein the first description comprises or is derived froman error function.
 67. The method according to claim 31, wherein thefirst description is a one dimensional description of temperaturedistribution about a node or element comprising an error function. 68.The method according to claim 31, wherein the solving step (d) utilizingthe conservation of mass and conservation of momentum equationscomprises the substeps of: (i) solving for pressure over at least partof the solution domain; and (ii) calculating velocity over at least partof the solution domain.
 69. The method according to claim 31, whereinthe solving step (d) comprises the substep of calculating temperaturebased on a convective heat transfer contribution, a conductive heattransfer contribution, and a viscous dissipation contribution.
 70. Themethod according to claim 41, wherein the solving step (e) utilizing theconservation of mass and conservation of momentum equations comprisesthe substeps of: (i) solving for pressure over at least part of thesolution domain; and (ii) calculating velocity over at least part of thesolution domain.
 71. The method according to claim 41, wherein thesolving step (e) comprises the substep of calculating temperature basedon a convective heat transfer contribution, a conductive heat transfercontribution, and a viscous dissipation contribution.
 72. The methodaccording to claim 26, wherein at least one of steps (d) and (e)comprises solving for fluidity over at least part of the solutiondomain.
 73. The method according to claim 26, wherein at least one ofsteps (d) and (e) comprises using a formulation derived from theconservation of mass and conservation of momentum equations.
 74. Themethod according to claim 26, wherein at least one of steps (d) and (e)comprises using at least one of a Nakano formulation, a Stokesformulation, and a Navier-Stokes formulation to solve for fluidity overat least part of the solution domain.
 75. The method according to claim31, wherein step (d) comprises solving for fluidity over at least partof the solution domain.
 76. The method according to claim 31, whereinstep (d) comprises using a formulation derived from the conservation ofmass and conservation of momentum equations.
 77. The method according toclaim 31, wherein step (d) comprises using at least one of a Nakanoformulation, a Stokes formulation, and a Navier-Stokes formulation tosolve for fluidity over at least part of the solution domain.
 78. Anapparatus for modeling injection of a fluid into a mold defining a threedimensional cavity, the apparatus comprising: a memory for storing codethat defines a set of instructions; and a processor for executing saidset of instructions to: (a) discretize a solution domain based on athree dimensional computer model defining a cavity; (b) at least one of:(i) solve for filling phase process variables over at least part of thesolution domain to provide respective filling phase solutions therefore;and (ii) solve for packing phase process variables over at least part ofthe solution domain using respective states of the process variables attermination of filling to provide respective packing phase solutionstherefor, wherein at least one of (i) and (ii) comprises: using a firstdescription of a distribution of a process variable about each of aplurality of nodes or elements within the solution domain; and using asecond description of the distribution of the process variable in atleast part of the solution domain comprising the plurality of nodes orelements, the second description comprising conservation of mass,conservation of momentum, and conservation of energy equations.